# Closure of Bessel Functions of the first kind

I need to use the Bessel functions of the first kind to solve some initial value problem. For this I need the closure equation

$$\int_0^\infty J_m(au)J_m(bu)u\,\text{d}u = \frac{\delta(a-b)}{a} \quad\quad \text{for}\quad a, b, m \in \mathbb{R}\quad \wedge \quad a, b >0 ,$$ which can also be found on the Wolfram functions webpage (see also Arfken and Weber, p.696, Morse and Feshbach, Section 6.3). However Mathematica (12.3.1.0) does some weird things:

In[12]:= Refine[Integrate[u*BesselJ[1, b u]* BesselJ[1, a u], {u, 0, Infinity}], a > 0 && b > 0 && a != b]

Out[12]= ConditionalExpression[0, a > b]

In[13]:= Integrate[u*BesselJ[1, 2 u]* BesselJ[1,  3 u], {u, 0, Infinity}]

During evaluation of In[13]:= Integrate::idiv: Integral of u BesselJ[1,2 u] BesselJ[1,3 u] does not converge on {0,\[Infinity]}.


So according to the first expression AND the closure equation the second integral should evaluate to 0. However the integral diverges. Since the closure relation is on the wolfram functions page I assume Mathematica should be able to apply it.

Some colleagues tried evaluating the expression in Mathematica 12.0 and get the correct result.

Is there something wrong in my code?

Edit:

As @yarchik has pointed out, the two results do not necessarily contradict each other. Since the the Dirac delta is a generalized function we cannot assign a specific value to $$\delta(a-b) \quad \text{for} \quad a \neq b.$$ In particular the delta distribution is defined to act on test functions $$f$$, that is we have to consider integrals of the form $$\int \delta(a-b) f(b) db$$ to assign some sort of value to it.

Of course in initial value problems integrals of this type are usually encountered. In my case I actually have the integral $$\int_0^\infty dk k \int_0^\infty du u J_m(ku)J_m(k'u) = \int_0^\infty dk k \frac{\delta(k-k')}{k} = 1,$$ which is well defined.

The problem was, that I tried to solve the inner integral on its own and expected it to be zero for $$a\neq b$$ which is a wrong assumption.

• Welcome to Mathematica.SE! I hope you will become a regular contributor. To get started, 1) take the introductory tour now, 2) when you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge, 3) remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign, and 4) give help too, by answering questions in your areas of expertise. Jul 15 at 10:41
• As I remember it, the question was asked and commented at this forum. Unfortunaly,I can't instantly find a reference. The result of NIntegrate[u*BesselJ[1, u]*BesselJ[1, 2 u], {u, 0, Infinity}] which is -0.283769 contradicts the statement of the question. – Jul 15 at 11:05
• I had a look, but couldn't find a similar question. The closure statement of Bessel functions is a well known result and is often used initial value problems. So there is no doubt that it is true (See for example p. 696 in Arfken, Weber or Morse and Feshbach, Section 6.3). I also found an answer on here, where the closure relation was used to determine initial conditions. The integrand is however heavily oscillating, so I assume that is why NIntegrate does give a wrong result. Jul 15 at 11:34
• @Mysterioso: First, NIntegrate  produces no warning/error. Second, Normal[Series[u*BesselJ[1, 2 u]*BesselJ[1, 2 u], {u, Infinity, 1}]] produces Cos[\[Pi]/4 + 2 u]^2/\[Pi] - ( 3 Cos[\[Pi]/4 + 2 u] Sin[\[Pi]/4 + 2 u])/(8 \[Pi] u) and this implies the divergence. Third, the command of Maple evalf(Int(u*BesselJ(1, u)*BesselJ(1,2* u), u = 0 .. infinity)) results in Float(undefined). Jul 15 at 12:00
• @user64494 So you are saying the Closure equation is wrong? Jul 15 at 12:04

There is a numerical way to check this Mathematica

i[a_, m_] := NIntegrate[
Exp[-u/100] BesselJ[m, a u]*BesselJ[m, 1/2 u] , {u, 0, Infinity}]


where I set $$b=1/2$$. Notice that some regularization, as is typical for expressions yielding generalized functions, is required. Now we can do some plotting

Plot[{i[a, 1], i[a, 2], i[a, 3], i[a, 4], i[a, 5]}, {a, -1, 1}]


It reveals two sharp maxima, not just one as in your equation

• Thanks for showing the nice way to analyze such expressions. You are right, $a$ and $b$ represent wave numbers and should actually be bigger than 0. I added this to the question. Allthough this should not impact the second integral of my question, since $a$ and $b$ are both chosen to be positive there. Do you have an idea why Mathematica does not use the relation it gave when calculating the general integral? Jul 16 at 8:04
• @Mysterioso Why do you think that one is more correct than another? The first general answer is correct only in some generalized sense. Its like integrating $\int_{-\infty}^\infty e^{ikx}dx$. If you put in some $k$, let us say $k=1$ you get $\int_{-\infty}^\infty e^{i x}dx$ clearly diverging. But if you consider it as a function of $k$ inside some other integral, you can perform the integration in some sense using $\delta(k)$. Jul 16 at 8:23
• Yeah, you are right. I confused the Dirac delta distribution with the Kronecker delta. I do not have much experience with using generalized functions in Mathematica. Thank you for your answer! Jul 16 at 8:34