Gaussian integral $\int_{0}^{\infty} dx \, x \, e^{-a^2 x^2} \left( \sqrt{(x - c)^2 + b^2} - \sqrt{(x + c)^2 + b^2} \right)$ I've been having trouble evaluating the Gaussian integral of the form
$$\int_{0}^{\infty} dx \, x \, e^{-a^2 x^2} \left( \sqrt{(x - c)^2 + b^2} - \sqrt{(x + c)^2 + b^2} \right) \, ,$$
which can be rearranged into integrals of the form
$$\int_{0}^{\infty} dx \, x^k \sqrt{x^2 + b^2} \left( e^{-a^2 (x + c)^2} \pm e^{-a^2 (x - c)^2} \right) \, ,$$
for $a, \, b, \, c > 0$ and $k \in \mathbb{Z}$. Does anyone know how to evaluate these expressions analytically? Or know any good integral tables on Gaussians I should refer to?
Edit: While browsing Gradshteyn and Ryzhik, I came across notation for my difference of square roots in the section on Bessel functions, see $\S 6.52$ and e.g. notation $6.522$ with $\ell_1 = \frac{1}{2} \left( \sqrt{(b + c)^2 + a^2} - \sqrt{(b - c)^2 + a^2} \right)$. I didn't realise functions of $\ell_1$ showed up in other contexts, although here $\ell_1$ appears in the integrand and as a function of $x$. Based on this connection though, I wonder if the integral above may evaluate to a Bessel function or combination thereof. Of course, this integral is still very non-trivial, so any help would be much appreciated.
 A: As @rtem Alexandrov commented, a formal expansion around $x=0$ would perfectly work and would probably converge fast since
$$\sqrt{(x - c)^2 + b^2} - \sqrt{(x + c)^2 + b^2}=-2 c+\frac{b^2 c}{x^2}+\frac{b^2 c \left(4 c^2-3 b^2\right)}{4x^4}+O\left(\frac{1}{x^6}\right)$$
Expanded as series around $x=0$, we have $$\sqrt{(x - c)^2 + b^2}=\sum_{n=0}^\infty \alpha_n x^n$$
$$\alpha_0=\sqrt{b^2+c^2}\qquad \alpha_1=-\frac{c}{\sqrt{b^2+c^2}} \qquad\alpha_n=\frac{c (2 n-3) \alpha_{n-1}-(n-3) \alpha_{n-2}}{n \left(b^2+c^2\right)} $$
$$\sqrt{(x + c)^2 + b^2}=\sum_{n=0}^\infty \beta_n x^n$$
$$\beta_0=\sqrt{b^2+c^2}\qquad \beta_1=\frac{c}{\sqrt{b^2+c^2}} \qquad\beta_n=-\frac{c (2 n-3) \beta_{n-1}+(n-3) \beta_{n-2}}{n \left(b^2+c^2\right)} $$
So, making the problem more general
$$I_k=\int_{0}^{\infty}  x^k \, e^{-a^2 x^2} \left( \sqrt{(x - c)^2 + b^2} - \sqrt{(x + c)^2 + b^2} \right) \,dx$$ Since $\alpha_{2n}=\beta_{2n}$, then
$$\large\color{red}{I_k=\frac 1{2a^{k+2}}\sum_{n=0}^\infty \frac{\alpha_{2n+1} -\beta_{2n+1}}{a^{2n}}\Gamma \left(n+1+\frac{k}{2} \right)}$$
For illustration purposes, trying with only ten terms and $a=1$,$b=2$, $c=3$, $k=1$ we get
$$-\frac{21175391684262887968287 }{14394292089850546320512}\sqrt{\frac{\pi }{13}}=\color{red}{-0.723176}66$$ while numerical integration gives $\color{red}{-0.72317639}$.
