Calculate the expectation of Inverse Bessel Process Simply put, I want to calculate the integral: 
$$(2\pi t)^{-3/2}\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}\frac{e^{-\frac{(x-a)^2+(y-b)^2+(z-c)^2}{2t}}}{\sqrt{x^2+y^2+z^2}}dzdydx$$
This integral is used to calculate the expecation of Inverse Bessel Process starts from $(a,b,c)$, and when $(a,b,c)=(1,0,0)$ the calculation is shown in a paper about local martingales and filtration shrinkage. 
The idea in the paper is to change of variable of polar coordinates, but the technique is  useful when $b=c=0$, however for arbitrary $b$ and $c$, I try the same technique, which integrates over $z$ and $y$ first, and run into integral like  
change of variable
$$y=r\cos{\theta},z=r\sin{\theta}$$
$$\int_{0}^{2\pi}\int_{0}^{\infty}\frac{e^{-\frac{(r\cos{\theta}-b)^2+(r\sin{\theta}-c)^2}{2t}}}{\sqrt{x^2+r^2}}rdrd\theta$$
change of variable again
$$\nu=\sqrt{x^2+r^2}$$
$$\int_{0}^{2\pi}\int_{0}^{\infty}e^{-\frac{\nu^2-x^2-2(b\cos{\theta}+c\sin{\theta})\sqrt{\nu^2-x^2}+b^2+c^2}{2t}}d\nu d\theta$$
Then I'm stuck here. I can't think of a change of variable that would further simplify this integral.
Anybody can show me how to do the original integral? You can ignore my attempt and use a new method.
 A: $\newcommand{\+}{^{\dagger}}%
 \newcommand{\angles}[1]{\left\langle #1 \right\rangle}%
 \newcommand{\braces}[1]{\left\lbrace #1 \right\rbrace}%
 \newcommand{\bracks}[1]{\left\lbrack #1 \right\rbrack}%
 \newcommand{\dd}{{\rm d}}%
 \newcommand{\isdiv}{\,\left.\right\vert\,}%
 \newcommand{\ds}[1]{\displaystyle{#1}}%
 \newcommand{\equalby}[1]{{#1 \atop {= \atop \vphantom{\huge A}}}}%
 \newcommand{\expo}[1]{\,{\rm e}^{#1}\,}%
 \newcommand{\floor}[1]{\,\left\lfloor #1 \right\rfloor\,}%
 \newcommand{\ic}{{\rm i}}%
 \newcommand{\imp}{\Longrightarrow}%
 \newcommand{\ket}[1]{\left\vert #1\right\rangle}%
 \newcommand{\pars}[1]{\left( #1 \right)}%
 \newcommand{\partiald}[3][]{\frac{\partial^{#1} #2}{\partial #3^{#1}}}
 \newcommand{\pp}{{\cal P}}%
 \newcommand{\root}[2][]{\,\sqrt[#1]{\,#2\,}\,}%
 \newcommand{\sech}{\,{\rm sech}}%
 \newcommand{\sgn}{\,{\rm sgn}}%
 \newcommand{\totald}[3][]{\frac{{\rm d}^{#1} #2}{{\rm d} #3^{#1}}}
 \newcommand{\ul}[1]{\underline{#1}}%
 \newcommand{\verts}[1]{\left\vert #1 \right\vert}%
 \newcommand{\yy}{\Longleftrightarrow}$
Let $\ds{\vec{r} \equiv \pars{x,y,z}\,,\quad\vec{R} \equiv \pars{a,b,c}.\quad}$ The integral becomes:

\begin{align}
\pp&=\pars{2\pi t}^{-3/2}
\int{\expo{-\verts{\vec{r} - \vec{R}}^{2}/\pars{2t}} \over r}\,\dd^{3}\vec{r}\,,\quad
\pars{~\mbox{The integration is over}\ {\mathbb R}^{3}~}
\end{align}

\begin{align}
\pp&=\pars{2\pi t}^{-3/2}
\int_{0}^{2\pi}\dd\phi\int_{0}^{\infty}\dd r\,r^{2}\,
{1 \over r}\int_{0}^{\pi}\dd\theta\,\sin\pars{\theta}
\expo{-\bracks{r^{2} - 2rR\cos\pars{\theta} + R^{2}}/2t}
\\[3mm]&=
{t^{-3/2} \over \root{2\pi}}\expo{-R^{2}/2t}\int_{0}^{\infty}\dd r\,r
\expo{-r^{2}/2t}\int_{0}^{\pi}\dd\theta\,\sin\pars{\theta}
\expo{rR\cos\pars{\theta}/t}
\\[3mm]&=
{t^{-3/2} \over \root{2\pi}}\expo{-R^{2}/2t}\int_{0}^{\infty}\dd r\,r
\expo{-r^{2}/2t}\braces{\left.
{\expo{rR\cos\pars{\theta}/t} \over -Rr/t}
\right\vert_{0}^{\pi}}
\\[3mm]&=
{t^{-3/2} \over \root{2\pi}}\expo{-R^{2}/2t}\,{t \over R}\int_{0}^{\infty}\dd r\,
\expo{-r^{2}/2t}\pars{-\expo{-rR/t} + \expo{rR/t}}
\\[3mm]&=
-\,{1 \over \root{2\pi t}}\,{\expo{-R^{2}/2t} \over R}\sum_{\sigma = \pm}
\sigma\int_{0}^{\infty}\expo{-r^{2}/2t - \sigma rR/t}\,\dd r
\\[3mm]&=
-\,{1 \over \root{2\pi t}}\,{\expo{-R^{2}/2t} \over R}\sum_{\sigma = \pm}
\sigma\int_{0}^{\infty}\exp\pars{-{\bracks{r + \sigma R}^{2} - R^{2}\over 2t}}
\,\dd r
\\[3mm]&=
-\,{1 \over \root{2\pi t}}\,{1 \over R}\sum_{\sigma = \pm}
\sigma\int_{\sigma R}^{\infty}\exp\pars{-r^{2}\over 2t}\,\dd r
=
{1 \over \root{2\pi t}}\,{1 \over R}\sum_{\sigma = \pm}
\sigma\int^{\sigma R}_{0}\exp\pars{-r^{2}\over 2t}\,\dd r
\\[3mm]&=
{1 \over \root{2\pi t}}\,{1 \over R}\,\root{2t}\,{\root{\pi} \over 2}\sum_{\sigma = \pm}
\sigma\,{2 \over \root{\pi}}\int^{\sigma R/\root{2t}}_{0}\exp\pars{-r^{2}}\,\dd r
\\[3mm]&=
-{1 \over 2R}{\rm erf}\pars{-\,{R \over \root{2t}}}
+
{1 \over 2R}{\rm erf}\pars{{R \over \root{2t}}}
\end{align}

\begin{align}
&\color{#0000ff}{\large(2\pi t)^{-3/2}\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}\frac{e^{-\frac{(x-a)^2+(y-b)^2+(z-c)^2}{2t}}}{\sqrt{x^2+y^2+z^2}}dzdydx}
\\[3mm]&\color{#0000ff}{\large =
{1 \over \root{a^{2} + b^{2} + c^{2}}}\,
{\rm erf}\,\pars{\vphantom{\Huge A^{A^{A}}}\root{a^{2} + b^{2} + c^{2} \over 2t}}}
\end{align}

A: The change of variable $(x,y,z)\to(u+a,v+b,w+c)$ shows that this integral is $E_0[\|B_t-A\|^{-1}]$ where $A=(a,b,c)$ and the process $(B_t)$ is a standard 3-dimensional Brownian motion starting at the origin. 
The invariance of the distribution of $B_t$ by rotations around the origin shows that this depends only on $\|A\|=\sqrt{a^2+b^2+c^2}$, that is, $E_0[\|B_t-A\|^{-1}]=E_0[\|B_t-U\|^{-1}]$, where $U=(\|A\|,0,0)$. Since the value of $E_0[\|B_t-U\|^{-1}]$ is in the paper, you are done.
