# An integral involving three Bessel functions

I am looking for a closed form for the following integral

$$I = \int_0^\infty \mathrm{d} x \ x \ J_0(ax) J_0(bx) J_1(cx)$$

which can be thought of as a particular case of the more general integral

$$I(n_1,n_2,n_3) = \int_0^\infty \mathrm{d} x \ x \ J_{n_1}(ax) J_{n_2} (bx) J_{n_3} (cx)$$

i.e. $I=I(0,0,1)$. I am aware of the following related results:

• A closed form, analytic result for I(0,0,0), see for instance Ref. [1].
• A closed form, analytic result for I(0,1,1), also in Ref. [1].
• A closed form, analytic result for the general case $I(n_1,n_2,n_3)$ when $n_1 + n_2 + n_3 = 0$, in Ref. [2].
• A general expression in terms of the hypergeometric function $F_4$, which would include I(0,0,1) as a particular case, which however requires $c>a+b$. This is in Ref. [3].
• The analytic/numerical approach of Ref. [4]

Still I cannot find a result for $I=I(0,0,1)$ for every value of $a$,$b$,$c$. I have verified numerically that the integral exists finite, even when the condition $c>a+b$ is not met.

So far I have tried expressing $J_1(cx)$ using recursion relations to try and recover other known integrals, or integrating by parts one or more Bessel functions, or expressing the Bessel functions as a series. All these approaches seems to make the integral more complicated.

Similarly to this [5] other question, the motivation is coupling of momenta in Quantum Mechanics.

[1] G.N. Watson, "A Treatise on the Theory of Bessel Functions", (Cambridge University Press, Cambridge), 1966.

[2] A. D. Jackson and L. C. Maximon, "Integrals of Products of Bessel Functions", SIAM J. Math. Anal. 3, 446 (1971).

[3] W.N. Bailey, "Some Infinite Integrals Involving Bessel Functions", Proceedings of the London Mathematical Society 40, 37 (1936).

• This integral is discussed in mathoverflow.net/questions/295425/… in a very different context. – Alexandre Eremenko May 18 '18 at 2:34
• @AlexandreEremenko I think the integral that comes up in the mean motion problem is $\int_0^\infty J_0(ax)\,J_0(bx)\,J_1(cx)\,dx$ (which can be evaluated exactly using ordinary and inverse trigonometric functions) whereas zakk's integral has an extra $x$ in the integrand. However, I would love to know more about the relation between these various integrals (and why they both come up). – Gro-Tsen May 18 '18 at 12:07
• @AlexandreEremenko: thanks, that very interesting! – zakk May 18 '18 at 13:54
• @Gro-Tsen: the two integrals (with and without the extra $x$) are connected by integration by parts, provided that you can calculate $I(0,0,2)$. – zakk May 18 '18 at 13:54
• So unless I am mistaken, at least when $|a-b|\leq c\leq a+b$ (or something), the integral $\int_0^\infty J_0(ax)\,J_0(bx)\,J_1(cx)\,dx$ is equal to $\frac{1}{c}\left(1-\frac{1}{\pi}\arccos\frac{c^2-a^2-b^2}{2ab}\right)$ by Al Kashi's theorem (=law of cosines) and the expression found in Weyl's 1938 Mean Motion paper mentioned by @AlexandreEremenko. Do you think you can get something useful on $\int_0^\infty x\,J_0(ax)\,J_0(bx)\,J_1(cx)\,dx$ from this? – Gro-Tsen May 19 '18 at 11:57

Without loss of generality we can set $c=1$; Mathematica returns a closed form for $a=b$.

For $|a|<1/2$ the result is

$$I= \int_0^\infty \mathrm{d} x \ x \ J_0(ax) J_0(ax) J_1(x)=\frac{4K(\alpha)}{\pi^2\sqrt{1-4a^2}}\left[2E(\alpha)-K(\alpha)\right]\qquad\qquad(\ast)$$ where $\alpha=\tfrac{1}{2}-\tfrac{1}{2}\sqrt{1-4a^2}$, and $E$ and $K$ are elliptic integrals (as defined here and here).

For $|a|>1/2$ Mathematica returns a Meijer-G function, $$I=\frac{2}{\sqrt\pi} G_{3,3}^{2,1}\left(\frac{1}{4 a^2}\left| \begin{array}{c} 1,1,1 \\ \frac{1}{2},\frac{3}{2},\frac{1}{2} \\ \end{array}\right. \right),$$

which equals the real part of $(\ast)$, so I would conclude that the elliptic-integral expression can be used for all $a$.

I tried to check this numerically, comparing $$A=\int_0^\infty \mathrm{d} x \ x \ J_0(ax) J_0(ax) J_1(x)$$ with $$B={\rm Re}\,\left[\frac{4K(\alpha)}{\pi^2\sqrt{1-4a^2}}\left[2E(\alpha)-K(\alpha)\right]\right],\;\;\alpha=\tfrac{1}{2}-\tfrac{1}{2}\sqrt{1-4a^2}.$$ The numerical evaluation of $A$ has rapid oscillations for $a<1/2$ as well as for $a>1/2$, and only the envelope is described by $B$. I don't know whether or not this is an artefact of the numerics.
[Zakk tells me they are an artefact, so it seems the elliptic-integral expression stands: A=B for $a<1/2$ and $a>1/2$.]

blue: numerical evaluation of A; gold: plot of B, both as a function of $a$

• In think zakk already mentioned this in his post (when pointing to Ref. [3]). – Mateusz Kwaśnicki May 17 '18 at 9:35
• Many thanks for the reply! Actually I didn't know this particular result Carlo Beenakker has reported here. However, in my problem I cannot assume $a=b < 1/2$. I am thinking about it: maybe it will be useful in some low-transferred-momentum limit... – zakk May 17 '18 at 9:41
• I added the $a=b>1/2$ case as well. – Carlo Beenakker May 17 '18 at 9:42
• Are these results from some paper/book? I may want to have a look at it... – zakk May 17 '18 at 9:44
• Thanks for editing the answer. I checked the same numerically, and by tuning NIntegrate[] precision one can remove the oscillations, so I would say they are just artifacts. Regarding the full problem (with $a \neq b$) it seems like the only way to go is Mellin transform leading to a Mellin-Barnes representation, leading to an infinite series. Not a closed form, but hopefully slightly better than the original integral. Thanks again! – zakk May 17 '18 at 15:36

Gradshteyn and Ryzhik, 6.578.1, solves the integral (for c=1) in terms of Gamma factors and the Appell F4 $F4(1/2, 3/2; 1,1, a^2,b^2)$ when $a>0, b>0$, and $a+b<1.$ In certain cases the F4 can be reduced a product of Gauss hypergeometric 2F1's or, if half integer in the 3rd and 4th arguments, a single 2F1. The given arguments result in a miss of the product case by 1 unit. That 1 unit cannot be made up with recursion formulas. I think it unlikely that the F4 representation will simplify.

However, the F4 is a double series and with either a or b small enough, should converge rapidly. Alternatively, there is an integral relationship for the F4 which gets rid of the oscillations, i.e., $$F4(1/2,3/2;1,1,a^2,b^2) = (Gamma Factors) * \int_0^\infty x \,I_0(ax)\,I_0(bx)\,K_1(x) dx.$$

• Many thanks! I think this is equivalent to the fourth bullet point in my question, however your advice on numerical calculation of F4 are very much appreciated, moreso since it looks like there's so simpler closer form! – zakk May 18 '18 at 0:15
• Unfortunately it seems like, when these integrals are applied to Quantum Mechanics, the condition $a+b<c$ is quite unphysical. Its converse, $a+b>=c$ would imply that $a$, $b$, and $c$ for a triangle, which can be interpreted as momentum conservation. – zakk May 18 '18 at 0:19