Double integral with Hankel transform

Solution 1:

This won't answer your question, but may be of some help. Most of the information below comes from Chapter 12 of Ronald Bracewell's The Fourier Transform and Its Applications, Second Edition, Revised, 1986.

The 2-D Fourier Transform

$$ F(u,v) = \mathscr{F}\left\{f(x,y)\right\}=\int_{-\infty}^\infty \int_{-\infty}^\infty f(x,y)\space e^{-2\pi i (xv+vy)} \space dx dy$$

when there is circular symmetry, such that $$f(x,y) = f(r)$$ $$r^2 = x^2+y^2$$ and thus $$F(u,v)=F(q)$$ $$q^2 = u^2+v^2$$

can be expressed as a zeroth order Hankel Transform, with a strictly reciprocal transform

$$F(q) = \mathscr{H}_0\left\{f(r)\right\} = 2\pi \int_0^\infty f(r) \space J_0(2\pi qr)\space r\space dr$$

$$f(r) = \mathscr{H}_0\left\{F(q)\right\} = 2\pi \int_0^\infty F(q) \space J_0(2\pi qr)\space q\space dq$$

Bracewell provides the following theorems using the above particular Hankel Transform convention

$$\begin{align*}\mathscr{H}_0\left\{f(ar)\right\} &= a^{-2}F\left(\dfrac{q}{a}\right) \quad\mbox{Similarity}\\ \\ \mathscr{H}_0\left\{f(r)+g(r)\right\} &= F\left(q\right)+G\left(q\right) \quad\mbox{Addition}\\ \\ \mathscr{H}_0\left\{ \int_0^{2\pi} \int_0^\infty f(r')g(R)r' \space dr'\space d\theta \right\} &= F\left(q\right)G\left(q\right) \quad\mbox{2D Convolution of Circularly Symmetric Functions}\\ R^2 = r^2 + r'^2 -2rr'\cos{\theta} &\\ \\ \int_0^\infty |f(r)|^2 r \space dr &= \int_0^\infty |F(q)|^2 q \space dq \quad \mbox{Rayleigh}\\ \\ \int_0^\infty f(r)g^*(r) r \space dr &= \int_0^\infty F(q)G^*(q) q \space dq \quad \mbox{Power}\\ \\ \mathscr{H}_0\left\{\dfrac{d}{dr}\left[rf(r)\right]\right\} &= -\dfrac{d}{dq}\left[qF(q)\right] \quad \mbox{Derivative}\\ \\ \mathscr{H}_0\left\{\dfrac{d}{dr}f(r)\right\} &= -\dfrac{d}{dq}\left[q\mathscr{H}_0\left\{r^{-1}f(r)\right\}\right] \quad \mbox{Derivative}\\ \\ \mathscr{H}_0\left\{r\dfrac{d}{dr}f(r)\right\} &= -q^{-1}\dfrac{d}{dq}\left[q^2F(q)\right] \quad \mbox{Derivative}\\ \\ 2\pi \int_0^{\infty} f(r)r \space dr &= F(0) \quad \mbox{Definite Integral}\\ \\\end{align*}$$

There is no Shift Theorem, although you can search online for papers by Natalie Baddour that introduce a "generalized shift" and theorem for 2D Fourier Transforms with circular symmetry.

Bracewell provides the following Hankel Transform pairs that may be useful for deriving Hankel Transforms for your functions $f(x)$ and $g(y)$ above:

$$\begin{align*}\mathscr{H}_0\left\{\dfrac{1}{r^2+a^2}\right\} &= 2\pi K_0(2\pi aq) \\ \\ \mathscr{H}_0\left\{e^{-ar}\right\} &= \dfrac{2\pi a}{\left(4\pi^2q^2 + a^2\right)^\frac{3}{2}} \\ \\ \mathscr{H}_0\left\{\dfrac{e^{-ar}}{r}\right\} &= \dfrac{2\pi}{\left(4\pi^2q^2 + a^2\right)^\frac{1}{2}} \\ \\ \mathscr{H}_0\left\{e^{-\pi r^2}\right\} &= e^{-\pi q^2} \\ \\ \mathscr{H}_0\left\{r^2 e^{-\pi r^2}\right\} &= \left(\dfrac{1}{\pi} - q^2\right)e^{-\pi q^2} \\ \\ \end{align*}$$

Bracewell's book also touches on 3D Fourier Transforms in cylindrical coordinates, which may or may not be useful to you.

The following online book has tables of Hankel Transforms:

https://authors.library.caltech.edu/43489/7/Volume%202.pdf

But you may need to transform the convention used. The English Wikipedia page on the Hankel Transform has a description on how to transform the convention.

Good luck.

Update

Using the above relations, I was able to develop the following expression for $I$:

$$\begin{align*}I &= \left(\left[\dfrac{a}{\left(2\pi r^2+a^2\right)^\frac{3}{2}} \ast K_0\left(b\sqrt{2\pi}r\right)\right] \ast g\left(\sqrt{2\pi}r\right)\right) \biggr{|}_{r=0}\\ \\ &= \left(\left[\dfrac{a}{\left(y^2+a^2\right)^\frac{3}{2}} \ast K_0\left(by \right)\right] \ast g\left(y\right)\right) \biggr{|}_{y=0}\\ \end{align*}$$

Where $*$ denotes 2D convolution of circularly symmetric functions.

The good news is that it is not necessary to figure out the Hankel Transform of $g(y)$. The bad news is that you have two 2D convolutions to compute. At least you only need the value at the origin, so maybe there is some savings there. You can also convert the 2D convolutions from polar coordinates back to cartesian coordinates, if you think there is some simplification there, but I suspect there is not.

Solution 2:

$\textbf {Edition of 05.12.2018}$

$$\color{brown}{\textbf{1. Integral representation of the Bessel function}}$$

Let us consider another way using the integral representation of the Bessel function $(69)$

$$J_0(z)=\dfrac1\pi\int\limits_0^\pi \cos(z\sin\theta)\,\mathrm d\theta=\dfrac2\pi\int\limits_0^{\pi/2} \cos(z\sin\theta)\,\mathrm d\theta.\tag{1.1}$$

Substitution $\sin\theta = t$ transforms it to $$J_0(z) = \dfrac2\pi\int\limits_0^1\dfrac{\cos(zt)}{\sqrt{1-t^2}}\mathrm dt.\tag{1.2}$$ This leads to the formula $$I = \dfrac2\pi\int\limits_0^1\dfrac{C(t)}{\sqrt{1-t^2}}dt,\tag{1.3}$$ where $$C(t) = \int\limits_0^\infty \int\limits_0^\infty f(x)g(y)xy\cos(xyt)\,\mathrm{dx}\,\mathrm{dy}.\tag{1.4}$$ Easy to see that formulas $(1.3)-(1.4)$ allow to replace Bessel function to cosine one.

$$\color{brown}{\textbf{2. Testing on $1D$ Hankel transform}}$$

By the definition, $$F(q)=\mathcal H(f(x))= 2\pi\int\limits_0^\infty xf(x)J_0(2\pi qx)\,\mathrm dx.\tag{2.1}$$

Identity $(1.2)$ allows to write $(2.1)$ in the form of $$F(q)= 4\int\limits_0^1\int\limits_0^\infty \dfrac{xf(x)\cos(2\pi qxt)}{\sqrt{1-t^2}}\,\mathrm dx\,\mathrm dt = 4\int\limits_0^1 \dfrac{\Phi(2\pi qt)}{\sqrt{1-t^2}}\mathrm dt,\tag{2.2}$$ where $$\Phi(p) = \int\limits_0^\infty xf(x)\cos(px)\,\mathrm dx.\tag{2.3}$$ Let $$f(x) = \dfrac {e^{-sx}}x,$$ then $$F(q)=\mathcal H\left(\dfrac {e^{-sx}}x\right) = 2\pi \int\limits_0^\infty e^{-sx}J_0(2\pi qx)\mathrm dx = 2\pi\mathcal L(J_0(2\pi qx)) =\dfrac{2\pi}{\sqrt{s^2+4\pi^2q^2}}.\tag{2.4}$$ On the other hand, $$\Phi(p) = \int\limits_0^\infty e^{-sx}\cos(px)\,\mathrm dx = \mathcal L(\cos(px)) = \dfrac s{s^2+p^2}\tag{2.5}$$ and $$F(q) = 4\int\limits_0^1 \dfrac{\Phi(2\pi qt)}{\sqrt{1-t^2}}\mathrm dt = 4\int\limits_0^1 \dfrac{s}{s^2+4\pi^2q^2t^2}\dfrac{\mathrm dt}{\sqrt{1-t^2}} = \dfrac{2\pi}{\sqrt{s^2+4\pi^2q^2}}\tag{2.6}$$ (see also Wolfram Alpha).

Results $(2.4)$ and $(2.6)$ are the same.

Testing is successful.

$$\color{brown}{\textbf{3. Integration via the given function $f(x)$}}$$

For the given function $f(x),$ $$\int\limits_0^\infty\dfrac{xe^{-kx}}{x^2+b^2}\,\mathrm dx = \int\limits_0^\infty\dfrac{\dfrac xbe^{-kb\frac xb}}{\left(\dfrac xb\right)^2+1}\,\mathrm d\dfrac xb = \int\limits_0^\infty\dfrac{xe^{-kbx}}{x^2+1}\,\mathrm dx = F_1(kb),\tag{3.1}$$ where $$F_1(s)=\left(\dfrac\pi2 - \mathrm{Si}(s)\right)\sin(s) - \dfrac12\mathrm{Ci}(s)\cos(s))\tag{3.2}$$ (see also Wolfram Alpha).

Then $$\int\limits_0^\infty\dfrac{xe^{-ax}\cos(ytx)}{x^2+b^2}\,\mathrm dx = \dfrac{F_1(-ab+iybt)+F_1(-ab-iybt)}2,$$ $$C(t) = \int\limits_0^\infty yg(y)\dfrac{F_1(-ab+iybt)+F_2(-ab-iybt)}2\,\mathrm{dy}.\tag{3.3}$$ And the next step looks too complex.

$$\color{brown}{\textbf{4. Integration via the given function $g(y)$}}$$

For the given function $g(y),$ $$\int\limits_0^\infty yg(y)\cos(xty)\,\mathrm dy = \int\limits_0^\infty ye^{-y} L_l \left( \frac{2L+1}{l+L+1} y \right) L_L \left( \frac{2l+1}{l+L+1} y \right)\cos(xty)\,\mathrm dy = \dfrac{G(1-ixt)+G(1+ixt)}2,\tag{4.1}$$ where $$G(s) = \int\limits_0^\infty ye^{-sy} L_l(y+py) L_L(y-py)\,\mathrm dy,\tag{4.2}$$ or $$G(s) = -\dfrac{\mathrm d}{\mathrm ds}\mathcal L\bigl(L_l(y+py) L_L(y-py)\bigr),\tag{4.3}$$ where $$p = \frac{L-l}{l+L+1}\tag{4.4},$$

Is known formula for the convolution $$\int\limits_0^y L_m(t)L_n(y-t)\,\mathrm dt = \int\limits_0^y L_{m+n}(t)\,\mathrm dt = L_{m+n}(y)-L_{m+n+1}(y)\tag{4.5}$$ But I can't apply it for obtained kind of the integrals.

However, for the given values $l$ and $L$ the integral $(4.2)$ can be calculated and presented in the form of

$$G(s) = \dfrac {P(s)}{s^{L+l+1}},\tag{4.6}$$ wherein degree of the polynomial $P(s)$ is less than $L+l+1.$

Then $$C(t) = \dfrac12 \int\limits_0^\infty xf(x)\left(\dfrac{P(1-ixt)}{(1-ixt)^{L+l+1}} + \dfrac{P(1+ixt)}{(1+ixt)^{L+l+1}}\right)\,\mathrm dx,\tag{4.7}$$ wherein $$xf(x)\left(\dfrac{P(1-ixt)}{(1-ixt)^{L+l+1}} + \dfrac{P(1+ixt)}{(1+ixt)^{L+l+1}}\right)=\left(\dfrac{C_1(t)x+C_2(t)}{x^2+b^2}+\sum_{k=1}^{L+l+1}\dfrac{A_k(t)}{(1-ixt)^k}+\sum_{k=1}^{L+l+1}\dfrac{B_k(t)}{(1+ixt)^k}\right)e^{-ax}.$$ Let us calculate the integrals $$\int\limits_0^\infty\dfrac{e^{-kx}}{x^2+b^2}\,\mathrm dx = \dfrac1b \int\limits_0^\infty\dfrac{e^{-kb\frac xb}}{\left(\dfrac xb\right)^2+1}\,\mathrm d\dfrac xb = \dfrac1b \int\limits_0^\infty\dfrac{e^{-kbx}}{x^2+1}\,\mathrm dx = \dfrac1b F_0(kb),\tag{4.8}$$ where $$F_0(s)=\mathrm{Ci}(s)\sin(s) + \left(\dfrac\pi2-\mathrm{Si(s)}\right)\cos(s)\tag{4.9}$$ (see also Wolfram Alpha), $$\int\limits_0^\infty\dfrac{e^{-kx}}{1\pm ixt} = \dfrac 1t \int\limits_0^\infty\dfrac{e^{-\frac kt xt}}{1\pm ixt}\,\mathrm d(xt) = \dfrac 1t \int\limits_0^\infty\dfrac{e^{-\frac kt x}}{1\pm ix}\,\mathrm dx = \dfrac 1t \int\limits_0^\infty\dfrac{1\mp ix}{x^2+1}e^{-\frac kt x}\,\mathrm dx,$$ $$\int\limits_0^\infty\dfrac{e^{-kx}}{1\pm ixt} = \dfrac1t F_0\left(\dfrac kt\right) \mp \dfrac it F_1\left(\dfrac kt\right),\tag{4.10}$$ $$H_m(s)=\int\limits_0^\infty\dfrac{e^{-sx}}{(1\pm ix)^m}\,\mathrm dx = -\dfrac1{m-1} \int\limits_0^\infty e^{-sx}\,\mathrm d(1\pm ix)^{-m+1}$$ $$= -\dfrac1{m-1} e^{-sx}(1\pm ix)^{-m+1}\bigg|_0^\infty + \dfrac s{m-1}\int\limits_0^\infty\dfrac{e^{-sx}}{(1\pm ix)^{m-1}}\,\mathrm dx = \dfrac 1{m-1}\left(1 + sH_{m-1}(s)\right),$$ $$H_m(s)=\sum\limits_{n=0}^{m-2}\dfrac{(m-n-2)!}{(m-1)!}s^n+\dfrac1{(m-1)!}(F_0(s)\mp F_1(s)),\tag{4.11}$$ $$\int\limits_0^\infty\dfrac{e^{-kx}}{(1\pm ixt)^m}dx = \dfrac1t H_m\left(\dfrac kt\right).\tag{4.12}$$

Obtained formulas allow to calculate $C(t)$ in the closed form, after what remains using of $(1.3).$

Therefore, the main problem is to calculate $1$D integral which contains Ci and Si special functions.

These are the reasons why this approach can be more successful than Hankel transform.