I think I can show you another way to prove the following claim :
Claim :
$$\mathbb{E}\left(\frac{X}{\|X\|}\right)=\frac{\Gamma(\frac{n+1}{2})}{(2\sigma^2)^{1/2}\Gamma(\frac{n+2}{2})}{}_1F_1\bigg(\frac{1}{2};\frac{n+2}{2};-\frac{\|\mu\|^2}{2\sigma^2}\bigg)\mu$$
Proof :
Using River Li's answer, we have
$$\begin{align}
\mathbb{E}\left(\frac{X}{\|X\|}\right)
&= \mathbb{E}\left(\frac{X}{(\|X\|^2)^{1/2}}\right)
\\\\&=\mathbb{E}\left(X \cdot \frac{1}{\Gamma(1/2)}\int_0^\infty \mathrm{e}^{-z\|X\|^2}z^{-1/2}\,dz \right)
\\\\&= \frac{1}{\sqrt{\pi}}\int_0^\infty \mathbb{E}\left(X\mathrm{e}^{-z\|X\|^2}\right)z^{-1/2}\,dz
\\\\&= \mu \cdot \frac{1}{\sqrt{\pi}}\int_0^\infty \frac{e^{-\frac{\|\mu\|^2 z}{2\sigma^2 z + 1}} \cdot z^{-1/2}}{(2\sigma^2 z + 1)^{(n+2)/2}}dz\end{align}$$
Letting $t=\frac{1}{2\sigma^2 z + 1}$, we have $\frac{dt}{dz}=-2\sigma^2 t^2$, so we obtain
$$\mathbb{E}\left(\frac{X}{\|X\|}\right)=\mu\cdot \frac{e^{-\frac{\|\mu\|^2}{2\sigma^2}}}{(2\sigma^2)^{1/2} \sqrt{\pi}}\int_0^1 e^{\frac{\|\mu\|^2t}{2\sigma^2}}t^{\frac{n-1}2}(1-t)^{-1/2}dt$$
According to here, we can say that if $\text{Re}(b)\gt \text{Re}(a)\gt 0$, then we have
$$\int_0^1e^{zt}t^{a-1}(1-t)^{-a+b-1}dt=\Gamma(a)\Gamma(b-a)\ _{1}\bar{F}_1(a;b;z)$$
So, taking $a=\frac{n+1}{2},b=\frac{n+2}{2}$ and $z=\frac{\|\mu\|^2}{2\sigma^2}$, we have
$$\mathbb{E}\left(\frac{X}{\|X\|}\right)=\mu\cdot\frac{e^{-\frac{\|\mu\|^2}{2\sigma^2}}\Gamma(\frac{n+1}{2})}{(2\sigma^2)^{1/2}}{}_1\bar{F}_1\bigg(\frac{n+1}{2};\frac{n+2}{2};\frac{\|\mu\|^2}{2\sigma^2}\bigg)$$
Since $_{1}\bar{F}_{1}(a;b;z)=\frac{1}{\Gamma(b)}\ _{1}F_{1}(a;b;z)$, we have
$$\mathbb{E}\left(\frac{X}{\|X\|}\right)= \mu\cdot\frac{e^{-\frac{\|\mu\|^2}{2\sigma^2}}\Gamma(\frac{n+1}{2})}{(2\sigma^2)^{1/2}\Gamma(\frac{n+2}{2})}{}_1F_1\bigg(\frac{n+1}{2},\frac{n+2}{2};\frac{\|\mu\|^2}{2\sigma^2}\bigg)$$
Since $_{1}F_1(a;b;z)=e^z\ _{1}F_1(b-a;b;-z)$, we finally get
$$\mathbb{E}\left(\frac{X}{\|X\|}\right)=\frac{\Gamma(\frac{n+1}{2})}{(2\sigma^2)^{1/2}\Gamma(\frac{n+2}{2})}{}_1F_1\bigg(\frac{1}{2};\frac{n+2}{2};-\frac{\|\mu\|^2}{2\sigma^2}\bigg)\mu.\ \blacksquare$$
Added :
There is an error in step 2 since it is wrong that
$$\int_{1}^{-1} z \left(-e^{a z}\right) \left(1-z^2\right)^{\frac{n}{2}-\frac{1}{2}} dz = \frac{1}{2} \sqrt{\pi }\ a\ \Gamma\left(\frac{n+1}{2}\right) {}_0F_1\left(\frac{n}{2}+2;\frac{a^2}{4}\right)$$
According to WolframAlpha, for $(a,n)=(2,1)$, we have
$$LHS=\frac{e^4+3}{4e^2}\not=\frac{3\sqrt{\pi}}{4}\times\frac{e^4+3}{4e^2}=RHS$$
My conjecture is
$$\int_{1}^{-1} z \left(-e^{a z}\right) \left(1-z^2\right)^{\frac{n}{2}-\frac{1}{2}} dz = \frac{1}{2}\sqrt{\pi}\ a\ \frac{\Gamma\left(\frac{n+1}{2}\right)}{\color{red}{\Gamma\left(\frac{n+4}{2}\right)}}\ {}_0F_1\left(\frac{n}{2}+2;\frac{a^2}{4}\right)\tag1$$
According to WolframAlpha, for small $n$, $(1)$ holds.
However, I don't know how to prove that the conjecture is true.
Best Answer
The confluent hypergeometric functions are related to the generalized Laguerre polynomials: \begin{align} F(-n;x;t)&=\frac{\Gamma(n+1)\Gamma(x)}{\Gamma(x+n)}L_n^{(x-1)}(t) \end{align} so \begin{equation} H_{n,m}(x)=\frac{n!m!\Gamma^2(x)}{\Gamma(x+n)\Gamma(x+m)}\int_0^\infty t^{x-1}e^{-t}\ln t L_n^{(x-1)}(t)L_m^{(x-1)}(t)\,dt \end{equation} The orthogonality relation for the Laguerre polynomials reads \begin{equation} \int_0^\infty t^{x-1}e^{-t}L_n^{(x-1)}(t)L_m^{(x-1)}(t)\,dt=\frac{\Gamma(n+x)}{n!}\delta_{n,m} \end{equation} It can be differentiated with respect to $x$ to obtain \begin{equation} \int_0^\infty t^{x-1}e^{-t}\ln t L_n^{(x-1)}(t)L_m^{(x-1)}(t)\,dt+\int_0^\infty t^{x-1}e^{-t}\frac{d}{dx}\left[L_n^{(x-1)}(t)L_m^{(x-1)}(t)\right]\,dt=\frac{\Psi(n+x)\Gamma(n+x)}{n!}\delta_{n,m} \end{equation} From the differentiation relation \begin{equation} \frac{d}{dx}L_n^{(x-1)}(t)=\sum_{k=0}^{n-1}\frac{L_k^{(x-1)}(t)}{n-k} \end{equation} and recognizing the definition of $H_{n,m}(x)$, we have thus \begin{align} \frac{\Gamma(n+x)\Gamma(m+x)}{n!m!\Gamma^2(x)}H_{n,m}(x)&+\sum_{k=0}^{n-1}\frac{1}{n-k}\int_0^\infty t^{x-1}e^{-t}L_k^{(x-1)}(t)L_m^{(x-1)}(t)\,dt\\ &+\sum_{k=0}^{m-1}\frac{1}{m-k}\int_0^\infty t^{x-1}e^{-t}L_k^{(x-1)}(t)L_n^{(x-1)}(t)\,dt\\ &=\frac{\Psi(n+x)\Gamma(n+x)}{n!}\delta_{n,m} \end{align} using the orthogonality relation and supposing that $n> m$, only one term in the sums survives, while there is no if $n=m$: \begin{equation} \frac{\Gamma(n+x)\Gamma(m+x)}{n!m!\Gamma^2(x)}H_{n,m}+\frac{1}{n-m}\frac{\Gamma(m+x)}{m!} \left( 1-\delta_{n,m} \right)=\frac{\Psi(n+x)\Gamma(n+x)}{n!}\delta_{n,m} \end{equation} which is the proposed expression.