$\int_{-\infty}^\infty \frac{e^{-y^2/2}}{((y+y_0)^2+x_0^2)^r} dy$ I have to estimate the integral $$\int_{-\infty}^\infty \frac{e^{-y^2/2}}{((y+y_0)^2+x_0^2)^r} dy,$$ for $r\in \mathbb{R}^+$. I am a little amazed that Sage and Wolfram Alpha have nothing to say about it, and that Gradshteyn-Ryzhik doesn't seem to have anything on it either; it feels like a rather natural integral, the denominator being a distance function.
Of course I realize that, if I just expand $1/((y+y_0)^2+x_0^2)^r$ into a Taylor series around $y=-y_0$ and integrate term-by-term, I am not going to get something convergent; what I get is an asymptotic formula. But what is the right order of magnitude of the error term when the formula gets cut off at the $k$th term?
For instance: let $f(y)=1/((y+y_0)^2+(x_0^2))^r$ and write $$f(y) = f(y_0) + \frac{(y-y_0)^2}{2} O^*(\max_t f''(t)).$$ Then we get an error term of size $O(1/x_0^{2 r + 3})$. If we go up to a higher-order approximation, we obtain an error term of the form $O(1/x_0^{2(r+k)+1})$ for higher $k$. But can one also give an error term that depends on $y_0$ and not just on $x_0$, and thus is better when $x_0$ is large and $y_0$ is much larger still?
Note: also asked on MSE.
Note 2: just to be clear - the easy estimate (obtained by a cheap version of Laplace) is 
$$\frac{\sqrt{2 \pi}}{(x_0^2+y_0^2)^r} + O^*\left(\frac{2 r \cdot \sqrt{2 \pi}}{e^{3/2} x_0^{2 r + 2}}\right).$$ That is fine if $x_0 \asymp y_0$, but we would like to do better than that in general, particularly in the tricky case $y_0\ggg x_0\ggg 1$. Here $O^*(\text{X})$ means "something of absolute value at most $X$".
 A: I would have liked to say this in comments, but the tex integrals make it impractical. This is not an answer but a reformulation. We can change the integrand slightly: 
$$I(r)=I(r,x_0,y_0)= \int _{-\infty} ^{+\infty} e^{-y^2} \frac{dy}{((y+y_0)^2+x_0^2)^r}.$$ Multiplying by $\Gamma (r)$ we see that this becomes (after a change of variables) 
$$ \Gamma (r)I(r)= \sqrt{2\pi}\int _0^{\infty} \frac{du}{u} u^r \frac{1}{\sqrt{(1+u)}} e^{-(ux_0^2+y_0^2\frac{u}{1+u})}.$$ Hopefully, the latter integral can be estimated more easily.  
A: For $r$ an integer $\geq 2$ you could start from the closed form expression for $r=2$ and then obtain larger integer values of $r$ by differentiation with respect to $x_0$. The result for $r=2$ is
$$
I_2=\int_{-\infty}^\infty \frac{e^{-y^2/2}}{((y+y_0)^2+x_0^2)^2} dy=\frac{\sqrt{2 \pi }}{2x_0^2}+\frac{\pi e^{x_0^2/2}}{4x_0^3 e^{y_0^2/2}}(w+\bar{w}),
$$

with $z=x_0+iy_0$, $w=(1-x_0z)\text{erfc}(z/\sqrt{2})\exp(i x_0 y_0)$.

For the large-$x_0$ asymptotics I follow the suggestion of Venkataramana:
$$I_r=\int_{-\infty}^\infty \frac{e^{-y^2/2}}{((y+y_0)^2+x_0^2)^r} dy=\frac{\sqrt{2\pi}}{\Gamma(r)}\int_0^\infty\frac{u^{r-1}  }{\sqrt{2 u+1} }e^{-u \left(\frac{y_0^2}{2 u+1}+x_0^2\right)}\,du$$
Now for $x_0\rightarrow\infty$ the integrand contributes mainly near $u=0$, so I may replace $2u+1\mapsto 1$, and then find 
$$I_r\rightarrow\sqrt{2 \pi } \left(x_0^2+y_0^2\right)^{-r}.$$
Here are some numerical checks, blue is the exact integral, gold the large-$x_0$ form, all plotted as a function of $x_0$ for different values of $y_0$ and $r$.





top row from left to right: $y_0=1,r=2$; $y_0=1,r=4$
bottom row from left to right: $y_0=10,r=2$; $y_0=10,r=4$.

A: Mainly for purposes of comparison, 
let me flesh out what I called "a cheap version of Laplace". Write $\sigma = 2 r$.
Choose $\rho\in (0,1)$. Let $g(y) = 1/(x_0^2+y^2)^{\sigma/2}$.
  Then $$g''(y) = \frac{-\sigma}{(x^2+y^2)^{\frac{\sigma}{2} + 1}} +
  \frac{\sigma \left(\frac{\sigma}{2}+1\right) \cdot 2 y^2}{(x^2+y^2)^{\frac{\sigma}{2} + 2}}$$
  and so $|g''(y)|\leq \sigma (\sigma+1)/(x^2+y^2)^{\sigma/2+1}$. 
Let $I$ be the interval
 $\lbrack (1-\rho) y_0,
  (1+\rho) y_0\rbrack$. Then, for $y\in I$,
  $|g''(y)|\leq \sigma (\sigma+1)/((1-\rho) l_0)^{\sigma+2}$,
where $l_0 = \sqrt{x_0^2+y_0^2}$,
  and so
  $$g(y) = g(y_0) + g'(y_0) (y- y_0) + O^*\left(c_0 (y-y_0)^2\right),$$
  where $c_0 = \sigma (\sigma+1)/2 ((1-\rho) l_0)^{\sigma+2}$.
  Thus, by cancellation,
  $$\int_I g(y) e^{-(y-y_0)^2/2} dy = 
  \int_I g(y_0) e^{-(y-y_0)^2/2} dy +
  O^*\left(\int_I c_0 (y-y_0)^2 e^{-(y-y_0)^2/2} dy\right).$$
Since $g'(y)<0$ for $y\geq 0$,
we also know that $g(y)<g(y_0) + c_0 (y-y_0)^2$ 
for $y>(1+\rho) y_0$. We conclude that
  $$\begin{aligned}\int_{(1-\rho) y_0}^{\infty} g(y) e^{-(y-y_0)^2/2} dy 
  &\leq g(y_0) \int_{-\infty}^\infty e^{-y^2/2} dy +
  c_0 \int_{-\infty}^\infty y^2 e^{-y^2/2} dy\\
  &= g(y_0) \sqrt{2\pi} + c_0 \sqrt{2\pi} 
  =  \left(1 + \frac{\sigma (\sigma+1)}{2 (1-\rho)^{\sigma+2} l_0^2}\right)
  \frac{\sqrt{2\pi}}{l_0^\sigma} .\end{aligned}$$
It remains to consider $y\leq (1-\rho) y_0$. Since $g(y)\leq g(0) = 1/x_0^\sigma$,
  $$\begin{aligned}\int_{-\infty}^{(1-\rho) y_0} g(y) e^{-(y-y_0)^2/2} dy &=
  \frac{1}{x_0^\sigma} \int_{-\infty}^{-\rho y_0} e^{-y^2/2} dy \\ &\leq
  \frac{1}{x_0^\sigma (\rho y_0)} \int_{-\infty}^{-\rho y_0} y e^{-y^2/2} dy =
  \frac{e^{-\rho^2 y_0^2/2}}{\rho x_0^\sigma y_0}.
  \end{aligned}$$
Thus we obtain
  $$\int_{-\infty}^\infty \frac{e^{-(y-y_0)^2/2}}{(x_0^2+y^2)^{\sigma/2}} dy
  \leq \left(1 + \frac{\sigma (\sigma+1)}{2 (1-\rho)^{\sigma+2} l_0^2}\right)
  \frac{\sqrt{2\pi}}{l_0^\sigma} +
  \frac{e^{-\rho^2 y_0^2/2}}{\rho x_0^\sigma y_0}$$
for any $0<\rho<1$.
Hardly very powerful or elegant, but I wonder: (a) is the above qualitatively optimal? That is, are the lesser-order terms of the right order? (b) can one give an even quicker proof of the same or a stronger bound?
