# 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?

• This is an L^2 inner product of square-integrable functions, so it's the same if you consider the integral of the product of the Fourier transforms of the two functions. FT of Gaussian = Gaussian, I expect a quick residue computation to give the FT of the other part as some sort of piecewise decaying exponential (decaying on either side of the piecewise thing). This might make it more manageable. – Marcus Aurelius May 9 at 6:00
• That's interesting. The FT of $1/(y^2+x_0^2)^r$ is basically a modified Bessel function $K_{r-1/2}(x_0 |t|)$ (times $(t/x_0)^{r-1/2}$, times a factor depending only on $r$). The inner product seems to get muddled, however, in that the displacement becomes, naturally, a phase, and thus it is actually more work to get a gain from it than in the original problem. – H A Helfgott May 9 at 7:49

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$$.

Result

I have derived the following quickly converging series expansion for the integral

$$f = \sum_{k=0}^\infty a(k,x_0,y_0,r)\tag{1a}$$

where

$$a(k,x_0,y_0,r) \\=(\frac{\sqrt{\pi}}{k! \Gamma(r)}) 2^{\frac{1}{2}-k} e^{-\frac{y_0^2}{2}} \left(y_0^2\right)^k \left(2^{-r} \text{B}(r,k-r+\frac{1}{2}) \, _1F_1\left(r;-k+r+\frac{1}{2};\frac{x_0^2}{2}\right)+2^{-k-\frac{1}{2}} x_0^{2 k-2 r+1} \Gamma \left(-k+r-\frac{1}{2}\right) \, _1F_1\left(k+\frac{1}{2};k-r+\frac{3}{2};\frac{x_0^2}{2}\right)\right)\tag{1b}$$

where $$B(x,y) =\frac{\Gamma(x) \Gamma(y)}{\Gamma(x+y)}$$ is the Euler beta function and $$_1 F _1$$ is a hypergeometric function.

In order to avoid typing work here is the Mathematica expression

1/(k! Gamma[r]) 2^(
1/2 - k) E^(-(y0^2/
2)) Sqrt[\[Pi]] (y0^2)^k (2^(-(1/2) - k) x0^(1 + 2 k - 2 r)
Gamma[-(1/2) - k + r] Hypergeometric1F1[1/2 + k, 3/2 + k - r,
x0^2/2] +
1/Gamma[1/2 + k] 2^-r Gamma[1/2 + k - r] Gamma[
r] Hypergeometric1F1[r, 1/2 - k + r, x0^2/2])


Numerical example:

For $${r = 1, x_0 = 1, y_0 = 1}$$ the first 5 terms are
$$\{0.99686, 0.261743, 0.0415358, 0.00495025, 0.000477208\}$$

and their sum is $$1.30557$$. The numerical integration gives $$fn \simeq 1.30561$$.

Derivation

Here's a brief sketch of the derivation.

We start from the formula

$$z^{-r} = \frac{1}{\Gamma (r)} \int_0^{\infty } t^{r-1} \exp (-t z) \, dt\tag{2}$$

Letting $$z = x_0^2+(y+y_0)^2$$ and interchanging the order of integration the $$y$$-integral can be done exactly to give

$$\int_{-\infty }^{\infty } \exp \left(-t \left(x_0^2+(y+ y_0)^2\right)-\frac{y^2}{2}\right) \, dy\\ = \frac{\sqrt{2 \pi } e^{-t \left(\frac{y_0^2}{2 t+1}+x_0^2\right)}}{\sqrt{2 t+1}}$$

Now the $$t$$-integral

$$\frac{1}{\Gamma (r)} \left(\sqrt{2 \pi } e^{-\frac{y_0^2}{2}}\right) \int_0^{\infty } t^{r-1} e^{-t x_0^2}\frac{ exp({\frac{y_0^2/2}{2 t+1}})}{\sqrt{2 t+1}} \, dt$$

can't be done exactly but we can expand the exponential function containing $$y_0^2 \frac{1}{2t+1}$$ in a power series, and then the coefficients can be integrated to give $$(1)$$.

Special cases

The case $$r=0$$

gives $$f(0,x_0,y_0) = \sqrt{2 \pi}$$

The case $$y_0 = 0$$

$$f(r,x_0,y_0=0) = 2^{\frac{1}{2}-r} \Gamma \left(\frac{1}{2}-r\right) \, _1F_1\left(r;r+\frac{1}{2};\frac{x_0^2}{2}\right)+\left(x_0^2\right)^{\frac{1}{2}-r} B\left(\frac{1}{2},r-\frac{1}{2}\right) \, _1F_1\left(\frac{1}{2};\frac{3}{2}-r;\frac{x_0^2}{2}\right)$$

Asmptotic behaviour

$$f(r,x_0 \simeq 0,y_0=0) \simeq \left(x_0^2\right)^{\frac{1}{2}-r} B\left(\frac{1}{2},r-\frac{1}{2}\right)+2^{\frac{1}{2}-r} \Gamma \left(\frac{1}{2}-r\right) +O\left(x_0^2\right)$$

$$f(r,x_0 \simeq \infty,y_0=0) \simeq \sqrt{2\pi} x_0^{-2r}$$

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) 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?

• In particular, I would like to see that $\rho$ disappear. – H A Helfgott May 9 at 11:05