Limit of a double integral

What is the $\varepsilon\to 0$ limit of the following double integral $$\int\limits_{-1}^1d\tau\;\sqrt{1-\tau^2}\;\tau\int\limits_0^\infty dq\;q^2e^{iq(\tau+i\varepsilon)}\;?$$ I was asked about this integral by my friend who got it in a physics research project. In fact this is just $n=1$ case of a more general integral $$\int\limits_{-1}^1d\tau\;(1-\tau^2)^{n/2}\;\tau\int\limits_0^\infty dq\;q^{n+1}e^{iq(\tau+i\varepsilon)}.$$ For $n$ even, the integral can be calculated rather simply by using $$\int\limits_0^\infty dq \, e^{iq(\tau + i\varepsilon)} \,=\, \frac{\varepsilon}{\tau^2 \,+\, \varepsilon^2} \,+\, \frac{i\tau} {\tau^2 \,+\, \varepsilon^2} \,\to\, \pi\delta(\tau) \,+\, \mathcal{P}\frac{i}{\tau},$$ because in this case only the delta-function contributes (or more precisely its derivatives, as far as the original integral is concerned). But for $n$ odd, only the principal value part contributes and things become messy.

                          Note added


I believe Igor Khavkine's method of using Hadamard regularization is indeed most suitable for this integral. By the way, an useful reference about Hadamard principal value is http://cms.math.ca/10.4153/CJM-1957-015-1 (A generalization of the Cauchy principal value, by Charles Fox, Canad. J. Math. 9(1957), 110-117). However, brute force $\epsilon$-regularization method also works, as described below.

We can transform the initial integral in the following way $$I_n=\int\limits_{-1}^1d\tau\;(1-\tau^2)^{n/2}\tau\int\limits_0 ^\infty q^{n+1}e^{iq(\tau+i\epsilon)}dq=$$ $$(-1)^{n+1}\frac{d^{n+1}} {d\epsilon^{n+1}}\int\limits_{-1}^1d\tau\;(1-\tau^2)^{n/2}\tau\int\limits_0 ^\infty e^{iq(\tau+i\epsilon)}dq=$$ $$i(-1)^{n+1}\frac{d^{n+1}}{d\epsilon^{n+1}}\int\limits_{-1}^1d\tau\; (1-\tau^2)^{n/2}\frac{\tau^2+\epsilon^2-\epsilon^2}{\tau^2+\epsilon^2}=$$ $$i(-1)^n\frac{d^{n+1}}{d\epsilon^{n+1}}\left[\epsilon^2\int\limits_{-1}^1 \frac{(1-\tau^2)^{n/2}}{\tau^2+\epsilon^2}\;d\tau\right].$$ Let us make the substitution $\tau=\sin{\theta}$, then $$I_n=i(-1)^n\frac{d^{n+1}}{d\epsilon^{n+1}}\left[ \epsilon^2\int\limits_{-\pi/2}^{\pi/2}\frac{\cos^{n+1}{\theta}} {\sin^2{\theta}+\epsilon^2}\;d\theta\right]\equiv i(-1)^n\frac{d^{n+1}}{d\epsilon^{n+1}}\left[\epsilon^2J_{n+1}(\epsilon) \right]. \;\;\;\; (1)$$ We have $$J_n(\epsilon)=\int\limits_{-\pi/2}^{\pi/2}\frac{\cos^n{\theta}} {\sin^2{\theta}+\epsilon^2}\;d\theta=\int\limits_{-\pi/2}^{\pi/2} \frac{\cos^{n-2}{\theta}(1-\sin^2{\theta}-\epsilon^2+\epsilon^2)} {\sin^2{\theta}+\epsilon^2}\;d\theta,$$ and hence $$J_n(\epsilon)=(1+\epsilon^2)J_{n-2}(\epsilon)-K_{n-2},\quad \quad \quad \quad (2)$$ where $$K_n=\int\limits_{-\pi/2}^{\pi/2}\cos^n{\theta}\;d\theta.$$ Using (2), we can prove by induction $$J_{2m}=(1+\epsilon^2)^mJ_0-\sum\limits_{i=0}^{m-1}(1+\epsilon^2)^i K_{2(m-1)-2i}, \quad \quad \quad \quad (3)$$ $$J_{2m+1}=(1+\epsilon^2)^mJ_1-\sum\limits_{i=0}^{m-1}(1+\epsilon^2)^i K_{2m-1-2i}. \quad \quad \quad \quad (4)$$ Therefore, if $n=2m-1$ is odd, we get from (1) and (3) $$I_{2m-1}=-i\frac{d^{2m}}{d\epsilon^{2m}}\left [\epsilon^2(1+\epsilon^2)^mJ_0(\epsilon)-\epsilon^2(1+\epsilon^2)^{m-1} K_0\right ], \quad \quad \quad \quad (5)$$ all other terms from (3) giving zero. But $$J_0(\epsilon)=\int\limits_{-\pi/2}^{\pi/2}\frac{d\theta}{\sin^2{\theta}+ \epsilon^2}=\frac{\pi}{\epsilon\sqrt{1+\epsilon^2}},$$ and $K_0=\pi$. Therefore (5) reduces to $$I_{2m-1}=(2m)!i\pi-i\pi\frac{d^{2m}}{d\epsilon^{2m}}\left[\epsilon (1+\epsilon^2)^{m-1/2}\right ],$$ which is the same as $$I_{2m-1}=(2m)!i\pi-\frac{i\pi}{2m+1}\frac{d^{2m+1}}{d\epsilon^{2m+1}} (1+\epsilon^2)^{m+1/2}. \quad \quad \quad \quad (6)$$ In the $\epsilon\to 0$ limit, the second term in (6) gives zero, because the Taylor expansion of $(1+\epsilon^2)^{m+1/2}$ for small $\epsilon$ contains only even powers of $\epsilon$. Therefore, finally $$\lim_{\epsilon\to 0}I_{2m-1}=(2m)!i\pi. \quad \quad \quad \quad (7)$$ If $n=2m$ is even, we get similarly $$I_{2m}=i\frac{d^{2m+1}}{d\epsilon^{2m+1}}\left[\epsilon^2(1+\epsilon^2)^m J_1(\epsilon)\right],$$ because no term containing $K$-factors survives after taking $(2m+1)$-th derivative. Using $$\int\limits_{-\pi/2}^{\pi/2}\frac{\cos{\theta}}{\sin^2{\theta}+\epsilon^2} \;d\theta=\frac{2}{\epsilon}\arctan{\left(\frac{1}{\epsilon}\right)},$$ we get $$I_{2m}=2i\frac{d^{2m+1}}{d\epsilon^{2m+1}}\left [\epsilon(1+\epsilon^2)^m \arctan{\left(\frac{1}{\epsilon}\right)}\right ]. \quad \quad \quad \quad (8)$$ In the $\epsilon\to 0$ limit, since $\arctan{\left(1/\epsilon\right)}$ and its derivatives are not singular in this limit, the application of the Leibniz rule to (8) gives $$\lim_{\epsilon\to 0}I_{2m}=2i\lim_{\epsilon\to 0} \arctan{\left(\frac{1}{\epsilon}\right)}\frac{d^{2m+1}}{d\epsilon^{2m+1}} [\epsilon(1+\epsilon^2)^m]=i\pi(2m+1)! \quad \quad \quad \quad (9)$$ All other terms vanish, because, for any $i\ge 1$, either $$\lim_{\epsilon\to 0}\frac{d^{2i}}{d\epsilon^{2i}} [\epsilon(1+\epsilon^2)^m]=0,$$ since binomial expansion of $\epsilon(1+\epsilon^2)^m$ contains only odd powers of $\epsilon$, or $$\lim_{\epsilon\to 0}\frac{d^{2i}}{d\epsilon^{2i}} \arctan{\left(\frac{1}{\epsilon}\right)}=-\lim_{\epsilon\to 0}\frac{d^{2i-1}}{d\epsilon^{2i-1}}\frac{1}{1+\epsilon^2}=0.$$ Equations (7) and (9) can be unified in the final result $$\lim_{\epsilon\to 0}I_n=i\pi (n+1)!. \quad \quad \quad \quad (10)$$

-

I believe these integrals can be evaluated directly in the $\varepsilon=0$ limit by interpreting the result of the inner integral as a distribution. $$\int\limits_0^\infty dq \, q^{n+1} \, e^{iq(\tau + i\varepsilon)} \to (-i)^{n+1}\frac{\partial^{n+1}}{\partial\tau^{n+1}} \left(\pi\delta(\tau) + \mathcal{P}\frac{i}{\tau}\right).$$ The rest of the integrand is $\sqrt{1-\tau^2}\tau$, which is odd, so that for odd $n$ only the second term contributes. Evaluating the derivatives gives this term as $$i(-i)^{n+1}\frac{\partial^{n+1}}{\partial\tau^{n+1}}\mathcal{P}\frac{1}{\tau} = i(-i)^{n+1}\frac{\partial^{n+2}}{\partial\tau^{n+2}}\log|\tau| = -i^n (n+1)! \mathcal{P} \frac{1}{\tau^{n+2}},$$ where the $\mathcal{P}$ now stands for finite part rather than principal value. This distribution can be defined simply as the distributional derivative as in the middle expression or it can be defined using the Hadamard regularization, which consists of introducing an integration cutoff (say $\sigma$) at $\tau=0$ and taking $\sigma\to0$ after subtracting the unique polynomial in $1/\sigma$ that makes the limit finite (presuming no $\log|\sigma|$ terms are present).

The original integral is \begin{align*} I(n) &= -i^n(n+1)!\int_{-1}^1 d\tau \, \mathcal{P} \frac{(1-\tau^2)^{n/2}}{\tau^{n+1}} \\ &= -2^{-n}i^n(n+1)! \int_{-1}^1 dt \, \mathcal{P} \frac{1}{t^{n+1}} \frac{(1-t^2)^{n+1}}{(1+t^2)} \\ &= 2^{-n}i(-1)^{\frac{n+1}{2}}(n+1)! \, J(n), \end{align*} where I've made the substitution $\tau=\frac{2t}{1+t^2}$, which rendered the integrand rational.

Now, the Hadamard regularization can be applied to evaluate the $J(n)$ integral for any $n$: $J(1)=-2\pi$, $J(3)=8\pi$, $J(5)=-32\pi$, ... I'll extrapolate this to $J(n) = (-1)^{\frac{n+1}{2}} 2^n \pi$, which gives the final answer $$I(n) = i\pi (n+1)! \, .$$

Incidentally, to implement the Hadamard regularization, it suffices to apply the "student integration rule": find an anti-derivative and evaluate the difference at the end-points. With this trick, the extrapolation for $J(n)$ can be justified by computing the generating function for $J(n)$: \begin{align} j(z) &= \sum_{n=0}^\infty z^n J(n-1) = \int_{-1}^1 dt\, \mathcal{P} \frac{1}{1+t^2} \frac{1}{1-z (1-t^2)/t} \\ &= \frac{\pi}{2(1+4z^2)} - \frac{\log(-1)z}{4z^2+1} , \end{align} where the integral was evaluated using computer algebra. The annoyingly vague $\log(-1)$ term actually doesn't matter, since for odd values of $n$, we only need the $z$-even part of the generating function: $$J(-1) + z^2 J(1) + z^4 J(3) + \cdots = \frac{\pi}{2(1+4z^2)} = \sum_{n=-1,1,3,5,\ldots} z^{n+1} (-1)^{\frac{n+1}{2}} 2^n \pi .$$

-

Switch the order of integrations:

Mathematica gives $$\int_{-1}^1 \tau\sqrt{1-\tau^2}e^{iq(\tau+i\epsilon)}=i\pi e^{-\epsilon q}{J_2(q)\over q}.$$ Hence you want the limit of $$\int_0^\infty e^{-\epsilon q}J_2(q)q\,dq,$$ for which Mathematica finds 2.

-