Approximating $\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?
 A: If you want to approximate near $y=-y_0$ you should use what you have. But for farther values of $y$ i.e. $|y+y_0|\gg0$, you should see instead that
\begin{align}\frac1{[(y+y_0)^2+x_0^2]^r}&=\frac1{(y+y_0)^{2r}}\frac1{\left[1+\left(\frac{x_0}{y+y_0}\right)^2\right]^r}\\&=\sum_{n=0}^\infty\binom{-r}n\frac{x_0^{2n}}{(y+y_0)^{2r+2n}}\end{align}
which goes to $0$ as $|y|\to\infty$.
A: Mainly for purposes of comparison, 
let me give a bound using what amounts to a cheap version of Laplace's method. (Cross-posted on MathOverflow.)
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?
