# Integrating the product of a sine with two Bessel functions

I have a difficult integral again, which I would really prefer to be able to express analytically:

$$Y_{nlL}=\int_0^1 \sin \pi n r J_0( \gamma_l r )J_0(\gamma_L r) dr$$

Here $n=1,2,3,4,\dots$ and $\gamma_l$ is the $l$th zero of $J_0(x)$.

The integrand looks nice, but it's highly oscillatory for large values of the parameters, and since there are $3$ of them, precomputing it numerically for all the cases up to some $N_m,L_m$ takes very long.

I have found a lot of examples with a single Bessel function and the limits up to $\infty$ (see this question for example), but I don't see how they help in this case.

A series expansion of one of the functions is not the best option either, since the expression for the terms would remain complicated, if they even exist.

If the community thinks I should, I will provide some work of my own, but I would really prefer a hint in the correct direction without the need to do a lot of useless manipulations.

In the end, it would be great to get an expression in terms of a finite number of any special functions Mathematica can handle. Fast converging single series is an option too of course, but I don't see how it could work with an oscillatory integrand.

The motivation for this: my attempts to approximate a much more difficult double integral have led to this. Using numerical pre-computation of $Y_{nlL}$, I have achieved a good accuracy for the original integral, but what's more important: a very fast evaluation.

One way I could see is to use one of the integral representations for the Bessel function:

$$J_0(a)=\frac{2}{\pi} \int_0^1 \frac{\cos a u}{\sqrt{1-u^2}}du$$

This could allow me to easily take the integral over $r$, but I would be left with a double integral instead, with the integrand still oscillatory.

I suppose I could explicitly try Levin's quadrature for highly oscillatory integrals, but it needs integrand in the form:

$$f(x) e^{i g(x)}$$

So if I had to transform both the Bessel functions into this form, I would again be left with a double integral for each step of the quadrature.

• Unlikely to be solved analytically; see' Analysis of propagation...' G. Tyler, J. Optical Soc Am.vol. 7 #7 July 1990. There they analyze a 3 Bessel function integral (yours is too; consider $J_{1/2}$) for case where two of the arguments are different. Yours is even harder with 3 distinct arguments and a finite (though large) upper limit. I once did something numerically like this, and if I remember I calculated many of the signed areas between the zeros, and stored them in an array. After a while the zeros become somewhat regular and then you can apply series acceleration techniques. – skbmoore Sep 20 '18 at 2:38
• @skbmoore, thank you for the reference and the advice – Yuriy S Sep 20 '18 at 10:32