I am interested in finding a general formula for the following integral: $$\int_0^\infty J_\nu(x)^3dx,\tag1$$ where $J_\nu(x)$ is the Bessel function of the first kind: $$J_\nu(x)=\sum _{n=0}^\infty\frac{(-1)^n}{\Gamma(n+1)\Gamma(n+\nu+1)}\left(\frac x2\right)^{2n+\nu}.\tag2$$ Mathematica gives the following result: $$\int_0^\infty J_\nu(x)^3dx=\frac1{\pi\,\nu}{_3F_2}\left(\frac12,\frac12-\nu,\frac12+\nu;1-\frac\nu2,1+\frac\nu2;\frac14\right)+\frac{\Gamma\left(-\frac\nu2\right)\Gamma\left(\frac{1+3\nu}2\right)\cos\frac{\pi\,\nu}2}{2^{\nu+1}\ \pi^{3/2}\ \Gamma(\nu+1)}{_3F_2}\left(\frac{1-\nu}2,\frac{1+\nu}2,\frac{1+3\nu}2;1+\frac\nu2,1+\nu;\frac14\right),\tag3$$ which can be significantly simplified for odd and half-integer values of $\nu$. The results at those points allow to conjecture another, simpler general formula: $$\int_0^\infty J_\nu(x)^3dx\stackrel?=\frac{\Gamma\left(\frac16\right)\ \Gamma\left(\frac16+\frac\nu2\right)}{2^{5/3}\ 3^{1/2}\ \pi^{3/2}\ \Gamma\left(\frac56+\frac\nu2\right)},\tag4$$ which agrees with $(3)$ to a very high precision for many different values of $\nu$. It also has an advantage that it is defined for all $\nu>-\frac13$ whereas $(3)$ is undefined at every even $\nu$ and requires calculating a limit at those points.

Is it possible to prove the formula $(4)$?

Mathematica is able to evaluate $(1)$ for even values of $\nu$ in terms of the Meijer G-function. Plugging those expressions into $(4)$ we get another form of the conjecture: $$G_{3,3}^{2,1}\left(\frac14\left|\begin{array}{c}2a,1,2-2a\\\frac12,1-a,a\\\end{array} \right.\right)\stackrel?=\frac{\Gamma\left(\frac16\right)\ \Gamma\left(\frac23-a\right)}{2^{5/3}\ 3^{1/2}\ \pi\ \Gamma\left(\frac43-a\right)}.\tag5$$ Incidentally, the case $a=\frac12$ would positively resolve my another question.

  • $\begingroup$ Did Mathematica evaluated the integral directly? $\endgroup$ May 28, 2013 at 15:23
  • $\begingroup$ Did you use Mathematica 9 to evaluate these integrals? I tried wolfram alpha, but no answer. $\endgroup$ May 28, 2013 at 15:30
  • $\begingroup$ @MhenniBenghorbal Yes, Mathematica 9 can evaluate them: Integrate[BesselJ[\[Nu], x]^3, {x, 0, Infinity}, Assumptions -> Re[\[Nu]] > -1/3]. $\endgroup$ May 28, 2013 at 20:15
  • $\begingroup$ @VladimirReshetnikov Maybe it's also worth asking these questions on math overflow? $\endgroup$
    – Ovi
    Jul 5, 2013 at 17:22
  • $\begingroup$ nice question .........+1 $\endgroup$
    – TShiong
    Mar 4, 2023 at 21:21

3 Answers 3


Thank you for posting this question, I enjoyed trying to answer it.

Start with the expression that Mathematica gave you and replace each argument $\frac14$ of a hypergeometric function with $\frac z4$, because we will be taking limits. I will call the two hypergeometric functions $Q_1(z)$ and $Q_2(z)$. Each term can be brought to a closed form by using identity 16.6.2 from the DLMF.

Setting $a=\frac12$, $b=1-\frac\nu2$, we get $$ Q_1(z) = (1-z)^{-\frac12} F\left(\frac16, \frac36,\frac 56; 1-\frac\nu2, 1+\frac\nu2; \frac{-27 z}{4(1-z)^3} \right), $$ and setting $a=\frac{1+3\nu}{2}$, $b=1+\frac\nu2$, we get $$ Q_2(z) = (1-z)^{-\frac{1+3\nu}{2}} F\left(\frac{1+3\nu}6, \frac{3+3\nu}{6}, \frac{5+3\nu}{6}; 1+\frac\nu2, 1+\nu; \frac{-27z}{4(1-z)^3} \right). $$ (Note that there are 6 possible identities to try per function, one for each possible choice of $a$ and $b$ from the parameters, so it helps to do this on a computer.)

The reason this works is that now the point $z=1$ is a singular point of the hypergeometric functions on the right hand side, and Mathematica will succeed in finding the limits as $z\to1$. The expression for the whole integral that you have is $$ Q = \frac{2^{\frac43}\pi^{\frac12}}{\Gamma(-\frac16)\Gamma(\frac56+\frac\nu2)\Gamma(\frac56-\frac\nu2)\sin\frac{\nu\pi}2} \left( -1 + 3^{-\frac{3\nu}2}\cos\left(\frac{\nu\pi}{2}\right) \frac{\Gamma(\frac{1+3\nu}{2}) \Gamma(\frac56-\frac\nu2)}{\Gamma(\frac{1+\nu}2)\Gamma(\frac56+\frac\nu2)} \right). $$ Call the large expression in brackets $A$, and then write $$ A = -1 + B 3^{-\frac{3\nu}{2}}\cos\frac{\pi\nu}{2}, \qquad B = \frac{\Gamma(\frac{1+3\nu}{2}) \Gamma(\frac56-\frac\nu2)}{\Gamma(\frac{1+\nu}2)\Gamma(\frac56+\frac\nu2)}. $$

Now, Mathematica will not simplify $A$ or $B$ on its own, so it needs help. Set $x=\frac16+\frac\nu2$, and use the multiplication formula to get $$ \frac{\Gamma(\frac{1+3\nu}2)}{\Gamma(\frac{1+\nu}{2})\Gamma(\frac56+\frac\nu2)} = \frac{\Gamma(3x)}{\Gamma(x+\frac13)\Gamma(x+\frac23)} = \frac{\Gamma(x)}{2\pi} 3^{3x-1/2}. $$ After this, $A$ simplifies to $$ A = -1 + \frac{\Gamma(\frac16+\frac\nu2)\Gamma(\frac56-\frac\nu2)}{2\pi}\cos\frac{\pi\nu}{2} = -1 + \frac{\cos\frac{\pi\nu}{2}}{2\sin(\frac\pi6+\frac{\pi\nu}{2})}, $$ where I've also used the reflection formula for $\Gamma(z)\Gamma(1-z)$ to get rid of the gamma functions. Some further amount of manual trigonometry yields $$ A = -\frac{\sqrt{3}}{2}\frac{\sin\frac{\pi\nu}{2}}{\sin(\frac\pi6 + \frac{\nu\pi}{2})}. $$

Finally, write $$ \frac{1}{\sin(\frac\pi6+\frac{\pi\nu}{2})} = \frac{\Gamma(\frac16+\frac\nu2)\Gamma(\frac56-\frac\nu2)}{\pi}, $$ and substitute back. Lots of things cancel, and the answer is $$ Q = -\frac{3^{1/2}2^{1/3}}{\pi^{1/2}} \frac{\Gamma(\frac16+\frac\nu2)}{\Gamma(-\frac16)\Gamma(\frac56+\frac\nu2)}. $$

This closed form is equivalent to the one you gave through the use of $\Gamma(\frac16)\Gamma(-\frac16)=-12\pi$.

P.S. I would also like to note that the integral $$ I(\nu,c) = \int_0^\infty J_\nu(x)^2 J_\nu(c x)\,dx $$ and its general form $$ \int_0^\infty x^{\rho-1}J_\nu(a x) J_\mu(b x) J_\lambda(c x)\,dx $$ appear in Gradshteyn and Ryzhik, and you can find a paper "Some infinite integrals involving bessel functions, I and II" by W. N. Bailey, which evaluates this integral in terms of Appell functions, but only in the case $c>2$ ($|c|>|a|+|b|$), which is where the $F_4$ Appell function converges. DLMF 16.16.6 actually gives a way to write this integral as $$ \frac{\Gamma(\frac{1+3\nu}{2})c^{-1-2\nu}}{\Gamma(1+\nu)^2\Gamma(\frac{1-\nu}{2})} \,\,\,{}_2F_1\left( \frac{1+\nu}{2}, \frac{1+3\nu}{2}; 1+\nu; x \right)^2, \qquad x = \frac{1-\sqrt{1-4/c^2}}{2}, $$ but the issue is that this is only correct for $c>2$, and the rhs is complex for $c<2$. Appell function would only be defined by analytic continuation in this case anyway, and I didn't find anything useful about non-principal branches of Appell or hypergeometric functions.

For $c>2$, Mathematica also gives the following: $$ I(\nu,c) = \frac{\Gamma(\frac{1+3\nu}{2})c^{-1-2\nu}}{\Gamma(\frac{1-\nu}{2})\Gamma(1+\nu)^2} \,\,\,{}_3F_2\left( \frac{1+\nu}{2}, \frac{1}{2}+\nu, \frac{1+3\nu}{2}; 1+\nu, 1+2\nu; \frac{4}{c^2} \right), $$ but this is incorrect when $c<2$.

  • 1
    $\begingroup$ Great answer, thanks! I awarded an additional bounty worth +100 to it. $\endgroup$ Jul 9, 2013 at 2:51
  • $\begingroup$ @VladimirReshetnikov Thank you! $\endgroup$
    – Kirill
    Jul 9, 2013 at 3:34
  • $\begingroup$ Could you expand your Gradshteyn and Ryzhik citation a bit? I'm apparently not seeing it (in my 1992 printing). $\endgroup$ Apr 12, 2014 at 22:52
  • $\begingroup$ +1 for this great answer $\endgroup$
    – TShiong
    Mar 4, 2023 at 21:22

I tried a while ago to evaluate this integral using a generalization of Ramanujan's Master Theorem referred to as the method of brackets.

(I've previously evaluated integrals using this method here and here.)

I even had brief email exchange with Dr. Karen T. Kohl about this particular integral, someone who has authored and coauthored papers on this method.

I'm going to show what happens if you use this method to evaluate the case $\nu =0$ .

Let's first assume that $a>2$.

The Bessel function of the first kind of order zero has the series representation $$J_{0}(az) = \sum_{m=0}^{\infty} \frac{(-1)^{m}}{m!} \, \frac{1}{\Gamma(m+1)} \left(\frac{a^{2}z^{2}}{4} \right)^{m}. $$

And the square of the Bessel function of the first kind of order zero has the series representation $$J_{0}(z)^{2} = {_1F_2} \left(\frac{1}{2};1,1;-z^{2} \right) = \sum_{n=0}^{\infty} \frac{ (-1)^{n}}{n!} \, \frac{\Gamma \left(n+ \frac{1}{2} \right)}{\Gamma \left(\frac{1}{2} \right)} \frac{1}{\Gamma(n+1)^{2}} \, z ^{2n}.$$

Therefore, the integral $\int_{0}^{\infty} J_{0}(ax) J_{0}(x)^{2} \, \mathrm dx $ has the corresponding bracket series $$\sum_{m=0}^{\infty} \sum_{n=0}^{\infty} \phi_{m,n} \, \left(\frac{a}{2} \right)^{2m} \frac{1}{\Gamma(m+1)} \frac{\Gamma \left(n+ \frac{1}{2} \right)}{\Gamma \left(\frac{1}{2} \right)} \frac{1}{\Gamma(n+1)^{2}} \, \langle 2m+2n+1 \rangle. $$

Since there are two indices but only one bracket, there are two cases to look at.

If we let $m$ be a free parameter, then the bracket vanishes if $n= -m- \frac{1}{2} $,

and we have the contribution $$S_{1} = \frac{1}{2} \sum_{m=0}^{\infty} \phi_{m} \left(\frac{a}{2} \right)^{2m} \frac{1}{\Gamma(m+1)} \frac{\color{red}{\Gamma(-m)}}{\Gamma \left(\frac{1}{2} \right)} \frac{1}{\Gamma\left(-m + \frac{1}{2}\right)^{2}} \, \Gamma \left(m+ \frac{1}{2} \right). $$ However, since the terms of this series are not finite, we discard it.

If we let $n$ be a free parameter, then the bracket vanishes if $m = - n - \frac{1}{2},$ and we have the contribution $$ \begin{align} S_{2} &= \frac{1}{2} \sum_{n=0}^{\infty} \phi_{n} \, \left( \frac{a}{2} \right)^{-2m-1} \frac{1}{\Gamma \left(-n+ \frac{1}{2} \right)} \frac{\Gamma \left(n+ \frac{1}{2} \right)}{\Gamma \left(\frac{1}{2} \right)} \frac{1}{\Gamma(n+1)^{2}} \Gamma \left(n+ \frac{1}{2} \right) \\ &= \frac{1}{\pi a } \sum_{n=0}^{\infty} \frac{1}{n!} \frac{\Gamma \left(n+ \frac{1}{2} \right)^{3}}{\Gamma \left(\frac{1}{2} \right)} \frac{1}{\Gamma(n+1)^{2}} \, \left(\frac{4}{a^{2}} \right)^{m} \\ &= \frac{1}{a} \, {_3F_2} \left(\frac{1}{2}, \frac{1}{2}, \frac{1}{2}; 1,1; \frac{4}{a^{2}} \right). \end{align}$$

Therefore, for $a>2$, we have $$ \begin{align} \int_{0}^{\infty} J_{0}(ax) J_{0}(x)^{2} \, \mathrm dx = S_{2} &= \frac{1}{a} \, {_3F_2} \left(\frac{1}{2}, \frac{1}{2}, \frac{1}{2}; 1,1; \frac{4}{a^{2}} \right) \\ &\overset{(1)}{=} \frac{1}{a} \, {_2F_1} \left(\frac{1}{4}, \frac{1}{4}; 1 ; \frac{4}{a^{2}} \right)^{2}. \end{align}$$

But according to the Wolfram Functions Site, $$ {_2F_1} \left(\frac{1}{4}, \frac{1}{4}; 1; \frac{4}{a^{2}} \right) = \frac{2}{\pi } \, K \left( \frac{1-\sqrt{1- \frac{4}{a^{2}}}}{2} \right), $$ where $K(m)$ is the complete elliptical integral of the first kind with parameter $m=k^{2}$.

This identity provides the analytical continuation of $${_2F_1} \left(\frac{1}{4}, \frac{1}{4}; 1; \frac{4}{a^{2}} \right) $$ outside the unit circle.

It's then natural to wonder if this means that $$\int_{0}^{\infty} J_{0}(x)^{3} \, \mathrm dx \overset{?}{=} \Re \left( {_2F_1} \left(\frac{1}{4}, \frac{1}{4}; 1; 4\right)^{2} \right)= \frac{4}{\pi^{2}} \, \Re \left( K \left(e^{-i \pi/3} \right)^{2} \right). $$

A numerical approximation of the integral suggests that this is indeed true, which would mean that $$\Re \left( K \left(e^{-i \pi/3} \right)^{2} \right) = \frac{\pi^{2}}{4} \frac{\Gamma \left(\frac{1}{6} \right)^{2}}{2^{5/3} 3^{1/2} \pi^{3/2} \Gamma \left(\frac{5}{6} \right)}.$$

There is a proposal at the end of this paper to modify the rule that says that all divergent series should be discarded when using this method. I have yet to see see this discussed in other papers.

$(1)$ https://en.wikipedia.org/wiki/Clausen%27s_formula

  • $\begingroup$ Very nice (+1). $\endgroup$ Mar 4, 2023 at 20:50
  • $\begingroup$ @AaronHendrickson Thanks. $\endgroup$ Mar 4, 2023 at 21:43

See section 6 in http://arxiv.org/abs/1007.0667 for that type of triple-product integrals.


You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .