# An integral involving Bessel functions and trigonometric functions.

Let $$x>0$$, $$y>0$$, $$\alpha>0$$ and $$m$$ be an integer bigger or equal to one. Consider a following integral:

$$$${\mathcal J}_{\alpha,m}(x,y):= \int\limits_0^\pi J_\alpha( x \cdot \sin(\theta) )\cdot e^{\imath y \cdot \cos(\theta)} \cdot \left[ \sin(\theta) \right]^m d\theta$$$$

where $$J_\alpha()$$ is the Bessel function of the first kind. Now by using the power series expansion of the Bessel function in the integrand then by integrating the the series term by term and then by resuming the resulting series we have obtained the following result:

$$$${\mathcal J}_{0,1}(x,y):= 2 \cdot j_0 \left( \sqrt{x^2+y^2}\right)$$$$

In here $$j_0()$$ is the spherical Bessel function. The code below verifies the result numerically. We have:

{x, y} = RandomReal[{0, 2}, 2, WorkingPrecision -> 50]; M = 50; t =.;
NIntegrate[
BesselJ[0, x Sin[th]] Exp[I y Cos[th]] Abs[Sin[th]], {th, 0, Pi},
WorkingPrecision -> 15]
1/2 Take[Accumulate[
2^(4/2)  Table[ ((-(x^2/(2 y)))^m)/m! Sqrt[\[Pi]/2] y^(-(1/2))
BesselJ[1/2 + m, y], {m, 0, M}]], -5] // MatrixForm
(*Here we identify the terms in the sum as spherical Bessel functions \
and we used the generating function identity from \
https://en.wikipedia.org/wiki/Bessel_function .*)
2 D[1/y Cos[Sqrt[y^2 - 2 y (t - x^2/(2 y))]], t] /. t :> 0
(2 Sin[Sqrt[x^2 + y^2]])/Sqrt[x^2 + y^2]


Now my question would be what is the result like for other values of $$n,m$$.

• You don't need the absolute value signs on $\sin(\theta)^m$ since $\sin(\theta)^m > 0$ for $\theta \in (0, \pi)$ Oct 1, 2020 at 17:14
• Yes, sorry. It is just that I was that the limits were different before and in that case I did need the modulus. Now I changed it to make it consistent some other approach. Oct 1, 2020 at 17:17
• Did you mean "resumming the resulting series"? May 4, 2023 at 12:58

Not a general answer, but the derivation of explicit representations for special cases.

We can modify the expression by changing $$\theta\to\pi-\theta$$ \begin{align} {\mathcal J}_{\alpha,m}(x,y)&= \int_0^\pi J_\alpha( x \sin(\theta) ) e^{\imath y\cos(\theta)} \sin^m(\theta) \,d\theta\\ &= \left( \int_0^{\pi/2}+\int_{\pi/2}^\pi\right) J_\alpha( x \sin(\theta) ) e^{\imath y\cos(\theta)} \sin^m(\theta) \,d\theta\\ &= \int_0^{\pi/2} J_\alpha( x \sin(\theta) ) e^{\imath y\cos(\theta)} \sin^m(\theta) \,d\theta+\int_0^{\pi/2}J_\alpha( x \sin(\theta) ) e^{-\imath y\cos(\theta)} \sin^m(\theta) \,d\theta\\ &=2\int_0^{\pi/2} J_\alpha( x \sin(\theta) )\cos( y\cos(\theta)) \sin^m(\theta) \,d\theta \end{align} Using the Bessel representation $$$$J_{-1/2}\left( y\cos\theta \right)=\sqrt{\frac{2}{\pi}}\frac{\cos( y\cos(\theta))}{\sqrt{y\cos(\theta)}}$$$$ we can express $$$${\mathcal J}_{\alpha,m}(x,y)=\sqrt{2\pi y}\int_0^{\pi/2} J_\alpha( x \sin(\theta) )J_{-1/2}\left( y\cos\theta \right) \sin^m(\theta) \cos^{1/2}\theta \,d\theta$$$$ A similar integral is tabulated (G&R 6.683.2): $$$$\int_0^{\pi/2} J_\nu( z_1 \sin\theta )J_{\mu}\left( z_2\cos\theta \right) \sin^{\nu+1}(\theta) \cos^{\mu+1}\theta \,d\theta=\frac{z_1^\nu z_2^\mu J_{\nu+\mu+1}\left( \sqrt{z_1^2+z_2^2} \right)}{\sqrt{\left( z_1^2+z_2^2 \right)^{\nu+\mu+1}}}$$$$ when $$\Re\nu>-1,\Re\mu>-1$$. By choosing $$\nu=\alpha,\mu=-1/2,z_1=x,z_2=y$$, if $$m=\nu+1$$, we obtain $$$${\mathcal J}_{\alpha,\alpha+1}(x,y)=\sqrt{2\pi}\frac{x^\alpha J_{\alpha+1/2}\left( \sqrt{x^2+y^2} \right)}{\left( x^2+y^2 \right)^{\alpha/2+1/4}}$$$$ When $$\alpha=0$$, we find $${\mathcal J}_{0,1}(x,y)=2j_0\left( \sqrt{x^2+y^2} \right)$$ as expected.

Other results may be obtained from the recurrences relations for the Bessel function. For example, using $$$$J_{\alpha}(z)=\frac{2(\alpha+1)}{z}J_{\alpha+1}(z)-J_{\alpha+2}(z)$$$$ by taking $$z=x\sin\theta$$, it comes $$$${\mathcal J}_{\alpha,\alpha+3}(x,y)=\frac{2(\alpha+1)}{x}{\mathcal J}_{\alpha+1,\alpha+2}(x,y)-{\mathcal J}_{\alpha+2,\alpha+3}(x,y)$$$$ Both terms of the rhs have an explicit representation from the above expression.

• @ Paul Enta Thank you very much for this. But there is a typo in your formula for ${\mathcal J}_{\alpha,\alpha+1}(x,y)$. The exponent in the denominator should be half of what it is now, i.e. $1/2(\alpha+1/2)$. Oct 2, 2020 at 11:54
• You're welcome. Sorry for the mistake which is corrected now. Oct 2, 2020 at 12:07

In here we derive the result for $$m = 2 \theta + \alpha + 1$$ where $$\theta \in {\mathbb N}$$. We start from the definition and we expand the integrand in a series and integrate term by term. We have:

$$\begin{eqnarray} {\mathcal J}_{\alpha,m}(x,y) &=& \sum\limits_{n=0}^\infty \frac{(-1)^n}{n! (n+\alpha)!} \cdot \left( \frac{x}{2} \right)^{2 n+\alpha} \cdot \int\limits_0^\pi [\sin(\theta)]^{2 n+m+\alpha} \cdot e^{\imath y \cos(\theta)} d\theta \\ &=& 2^{(m-\alpha)/2} \sqrt{\pi} x^\alpha \sum\limits_{n=0}^\infty \frac{1}{n! (n+\alpha)!} \cdot (n + (\alpha+m-1)/2 )! (-\frac{x^2}{2 y} )^n \cdot \frac{J_{n+(\alpha+m)/2}(y)}{y^{(\alpha+m)/2}} \\ &=& 2^{(m-\alpha)/2} \sqrt{\pi} x^\alpha \sum\limits_{n=0}^\infty \frac{1}{n!} (n+\alpha+1)^{(\theta)} (-\frac{x^2}{2 y} )^n \cdot \frac{J_{n+(\alpha+m)/2}(y)}{y^{(\alpha+m)/2}} \\ &=& 2^{(m-\alpha)/2} \sqrt{\pi} x^\alpha \sum\limits_{n=0}^\infty \frac{1}{n!} (n+\alpha+1)^{(\theta)} (-\frac{x^2}{2 y} )^n \cdot \left.\frac{d^n}{d t^n} \frac{J_{(\alpha+m)/2}(\sqrt{y^2-2 y t})}{(\sqrt{y^2-2 y t})^{(\alpha+m)/2}} \right|_{t=0} \\ &=& 2^{(m-\alpha)/2} \sqrt{\pi} x^\alpha \left. \frac{d^\theta}{d \xi^\theta} \xi^{\alpha+\theta} \frac{J_{(\alpha+m)/2}(\sqrt{y^2 + \xi x^2})}{(\sqrt{y^2 + \xi x^2})^{(\alpha+m)/2}}\right|_{\xi=1} \\ &=& 2^{(m-\alpha)/2} \sqrt{\pi} x^\alpha \sum\limits_{q=0}^\theta \sum\limits_{p=q}^\theta \binom{\theta}{p} (\alpha+p+1)^{(\theta-p)} {\mathfrak C}_q^p \frac{x^{2 p}}{2^p \sqrt{x^2+y^2}^{(\alpha+m)/2+2p-q}} \cdot \left. \frac{d^q}{d u^q} J_{(\alpha+m)/2}(u) \right|_{u= \sqrt{x^2+y^2}} \end{eqnarray}$$

In the second line from the top we integrated over the angle by substituting for $$\cos(\theta)$$ and then expanding the exponential in a series integrating term by term and identifying the result as an expansion of a Bessel function. In the third line we introduced a Pochammer symbol $$n^{(\theta)} = \Gamma(n+\theta)/\Gamma(n)$$ and in the forth line we made use of a remarkable differential identity which is a generalization of equation 10.1.40 in http://people.math.sfu.ca/~cbm/aands/page_439.htm . Finally in the fifth line we re-summed the series owing to the fact that it is just an appropriate Taylor expansion about the origin. Finally in the sixth line we used the chain rule to evaluate the derivative in question and in doing so we defined coefficients $${\mathfrak C}^p_q := \sum\limits_{j=q}^p (-1)^{p-j} (2(p-j)-1)!! \binom{2p-1-j}{j-1} \binom{j}{q} (-(\alpha+m)/2)_{(j-q)}$$ for $$q=0,\cdots,\theta$$ and $$p=q,\cdots,\theta$$. The result can be further simplified using differential identities for Bessel function from here

 {x, y, alpha, m} =
RandomReal[{0, 2}, 4, WorkingPrecision -> 50]; M = 10; t =.;
theta = RandomInteger[{0, 10}];
m = 2 theta + alpha + 1;
myCs = Table[
Sum[(-1)^(p - j) (2 (p - j) - 1)!! Binomial[2 p - 1 - j,
2 (p - j)] Binomial[j,
q] Pochhammer[-((alpha + m)/2) - (j - q) + 1, j - q], {j, q,
p}], {q, 0, theta}, {p, q, theta}];

NIntegrate[
BesselJ[alpha, x Sin[th]] Exp[I y Cos[th]] Sin[th]^m, {th, 0, Pi},
WorkingPrecision -> 15]
2^(1/2 (-alpha + m)) Sqrt[\[Pi]] x^alpha Take[
Accumulate[
Table[(n + (alpha + m - 1)/2)!/(n! (n + alpha)!) (-(x^2/(2 y)))^
n BesselJ[1/2 (alpha + m) + n, y]/y^(1/2 (alpha + m)), {n, 0,
M}]], -5] // MatrixForm
2^(1/2 (-alpha + m)) Sqrt[\[Pi]] x^alpha Take[
Accumulate[
Table[1/n! Pochhammer[
n + alpha + 1, (m - alpha - 1)/2] (-(x^2/(2 y)))^
n (D[BesselJ[1/2 (alpha + m),
Sqrt[y^2 - 2 y t]]/(Sqrt[
y^2 - 2 y t])^(1/2 (alpha + m)), {t, n}] /. t :> 0), {n,
0, M}]], -5] // MatrixForm
2^(1/2 (-alpha + m)) Sqrt[\[Pi]] x^alpha D[
xi^(alpha + theta) BesselJ[1/2 (alpha + m),
Sqrt[y^2 + xi x^2]]/(Sqrt[y^2 + xi x^2])^(1/2 (alpha + m)), {xi,
theta}] /. xi :> 1
(*Sqrt[2 Pi] x^alpha \
BesselJ[alpha+1/2,Sqrt[x^2+y^2]]/(Sqrt[x^2+y^2])^(alpha+1/2)*)
2^(1/2 (-alpha + m)) Sqrt[\[Pi]] x^alpha Sum[
Binomial[theta, p] Pochhammer[alpha + p + 1,
theta - p] myCs[[1 + q, p - q + 1]] x^(2 p)/(
2^p (Sqrt[x^2 + y^2])^(((alpha + m)/2) + 2 p - q)) (
D[BesselJ[nu, u], {u, q}] /. {u :> Sqrt[x^2 + y^2],
nu :> ((alpha + m)/2)}), {q, 0, theta}, {p, q, theta}]


Update: The result above above can be generalized to arbitrary values of $$m \in {\mathbb R}$$ except now we are left with an infinite sum which however converges fast as we will demonstrate numerically below. But first let us state the result. We take sum number $$M \in {\mathbb N}$$ and $$M\ge 1$$ and and we have:

$$\begin{eqnarray} &&{\mathcal J}_{\alpha,m}(x,y) = 2^{\frac{(m-\alpha)}{2}} \sqrt{\pi} x^\alpha\\ && \sum\limits_{q=0}^M \sum\limits_{p=q}^M \binom{\frac{(m-\alpha-1)}{2}}{p} (\alpha+p+1)^{(\frac{(m-\alpha-1)}{2}-p)} {\mathfrak C}_q^p \frac{x^{2 p}}{2^p \sqrt{x^2+y^2}^{\frac{(\alpha+m)}{2}+2p-q}} \cdot \left. \frac{d^q}{d u^q} J_{(\alpha+m)/2}(u) \right|_{u= \sqrt{x^2+y^2}} \quad (ii) \end{eqnarray}$$

We claim that the series converges fast. Indeed we took $$M=20$$ and using the code below we sampled random values of $$x,y,\alpha,m \in (0,10)$$ twenty times and every time the numerical integral matched the series $$(ii)$$ at least to twelve decimal digits precision.

{x, y, alpha, m} =
RandomReal[{0, 10}, 4, WorkingPrecision -> 50]; M = 20;
NIntegrate[
BesselJ[alpha, x Sin[th]] Exp[I y Cos[th]] Sin[th]^m, {th, 0, Pi},
WorkingPrecision -> 15]
ll = 2^(1/2 (-alpha + m)) Sqrt[\[Pi]] x^alpha Table[
Binomial[(m - alpha - 1)/2, p] Pochhammer[
alpha + p + 1, (m - alpha - 1)/2 -
p] Sum[(-1)^(p - j) (2 (p - j) - 1)!! Binomial[2 p - 1 - j,
2 (p - j)] Binomial[j,
q] Pochhammer[-((alpha + m)/2) - (j - q) + 1, j - q], {j, q,
p}] x^(2 p)/(
2^p (Sqrt[x^2 + y^2])^(((alpha + m)/2) + 2 p - q)) (
D[BesselJ[nu, u], {u, q}] /. {u :> Sqrt[x^2 + y^2],
nu :> ((alpha + m)/2)}), {q, 0, M}, {p, q, M}];
N[Take[Accumulate[Flatten[ll]], -5], 15] // MatrixForm


Update 1: Now we will generalize the result. Let us take as before $$\alpha, x, y \in {\mathbb R}_+$$ and also take $$m_1,m_2 \in {\mathbb N}_+$$. We define a following integral:

$$\begin{eqnarray} {\mathcal J}_{\alpha,(m_1,m_2)} (x,y) := \int\limits_0^\pi J_\alpha(x \sin(\theta)) \cdot e^{\imath y \cos(\theta)} \cdot [\sin(\theta)]^{m_1} \cdot [\cos(\theta)]^{m_2} d\theta \end{eqnarray}$$

Now define certain constants . Firstly we have:

% nu <-- \frac{\alpha+m_1}{2} $$\begin{eqnarray} {\mathcal C}^{(l)}_{n,j}(\nu) := \sum\limits_{q=j+l}^n \sum\limits_{k=j+l}^n (-1)^{q+n+k+j-1} \binom{n}{q} \left( - \nu \right)^{(n-q)} \cdot \binom{k-j-l}{j} \cdot \frac{2^{2k-q} \cdot (2 k-q+1)^{(2q-2k)}}{(q-k)! \left( \nu +k-j\right)^{(2j+l-k)}} \end{eqnarray}$$

for $$l=0,1$$ and $$n=0,\cdots,M$$ and $$j=0,\cdots,n-l$$ . We note that the constant $${\mathcal C}^{(.)}_{.,.}(\nu)$$ above is a polynomial in the variable $$\nu$$.

Secondly we have: %nu <-- (\alpha+m_1)/2 $$\begin{eqnarray} {\mathfrak C}^p_q(\nu) := \sum\limits_{j=q}^p (-1)^{p-j} (2(p-j)-1)!! \binom{2p-1-j}{j-1} \binom{j}{q} (-\nu)_{(j-q)} \end{eqnarray}$$ for $$q=0,\cdots,M$$ and $$p=q,\cdots,M$$. Again we note that the constant $${\mathfrak C}^{.}_{.}(\nu)$$ is a polynomial in the variable $$\nu$$.

Now let us define the following quantities below. We have:

$$\begin{eqnarray} d_1(m,j)&:=& 2 \left(\left\lfloor \frac{m-3}{4}\right\rfloor +\left\lfloor \frac{m-2}{4}\right\rfloor +\left\lfloor \frac{m}{4}\right\rfloor +2\right) \cdot 1_{j=0} + \left(0 \vee 2 \left(\left\lfloor \frac{m-2}{4}\right\rfloor +\left\lfloor \frac{m-1}{4}\right\rfloor +\left\lfloor \frac{m}{4}\right\rfloor +1\right)\right) \cdot 1_{j=1} \\ d_2(m,j)&:=& \left( \left\lfloor \frac{m-1}{2}\right\rfloor +3 \left\lfloor \frac{m}{2}\right\rfloor +1 \right) \cdot 1_{j=0} + \left(0 \vee \left(3 \left\lfloor \frac{m-1}{2}\right\rfloor +\left\lfloor \frac{m}{2}\right\rfloor +1\right)\right) \cdot 1_{j=1} \end{eqnarray}$$

And then we define coefficients $$(A^{m}_{j,k_1,k_2} )_{m=0,j=0,k_1=0,k_2=0}^{M,1,d_1(m,j),d_2(m,j)}$$ recursively as follows:

$$\begin{eqnarray} &&A^{m+1}_{j,k_1,k_2}=\\ && \left(\alpha+j-2 m-1_{j=0} \nu + 1_{j=1}(\nu-1) \right) \cdot A^{m}_{j,k_1,k_2-1} +\\ && (k_2+1) \cdot A^{m}_{j,k_1-2,k_2+1} +\\ &&(k_2-1) \cdot A^{m}_{j,k_1,k_2-1} +\\ &&1_{j=1} \cdot A^{m}_{j-1,k_1,k_2-1} + \\ &&(-1_{j=0}) \cdot A^{m}_{j+1,k_1-2,k_2-1} + \\ &&(-1_{j=0}) \cdot A^{m}_{j+1,k_1,k_2-3} \end{eqnarray}$$

subject to $$A^{0}_{0,0,0}= 1$$.

$$\begin{eqnarray} &&{\mathcal J}_{\alpha,(m_1,m_2)} (x,y) =(-\imath)^{m_2} \cdot 2^{\frac{m_1-\alpha}{2}} \cdot \sqrt{\pi} \cdot x^\alpha \cdot \\ &&\sum\limits_{q=0}^M \sum\limits_{p=q}^M \binom{\frac{m_1-\alpha-1}{2}}{p} (\alpha+p+1)^{(\frac{m_1-\alpha-1}{2}-p)} \cdot \left. {\mathfrak C}^{(p)}_q(\nu) \right|_{\nu= \frac{\alpha+m_1}{2}} \cdot (\frac{x^2}{2})^p \cdot \\ && \sum\limits_{l=0}^1 \sum\limits_{j=0}^{q-l} \sum\limits_{j_1=0}^{1} \sum\limits_{k_1=0}^{\lfloor d_1(m_2,j_1)/2 \rfloor} \sum\limits_{k_2=0}^{\lfloor (m_2-j_1)/2 \rfloor} (-1)^{l-1} \cdot \left. {\mathcal C}^{(l)}_{q,j}(\nu) \right|_{\nu= \frac{\alpha+m_1}{2}} \cdot \frac{1}{2^{2j+l}} \cdot \left. A^{(m_2)}_{j_1,2 k_1,m_2+2 k_2-2 k_1}(A,\nu) \right|_{ \left( \begin{array}{rrr} A&=& (2 j+l -(\alpha+m_1)/2-2 p +j_1)/2 - m_2 \\ \nu&=& (\alpha+m_1)/2-l \end{array} \right) } \\ && \cdot x^{2 k_1} y^{m_2+2 k_2-2 k_1} \left( x^2+y^2\right)^{(2j +l - \frac{\alpha+m_1}{2} - 2 p+j_1)/2 - m_2} \cdot J_{\frac{\alpha+m_1}{2} - l - j_1}\left( \sqrt{x^2+y^2}\right) \end{eqnarray}$$

{x, y, alpha} = RandomReal[{0, 10}, 3, WorkingPrecision -> 50]; M = 20;
{m1, m2} = RandomInteger[{1, 3}, 2]; xi =.; m2 = 2 m2;
NIntegrate[
BesselJ[alpha, x Sin[th]] Exp[I y Cos[th]] Sin[th]^m1 Cos[th]^
m2, {th, 0, Pi}, WorkingPrecision -> 15]
CC = Table[
Sum[(-1)^(p - j) (2 (p - j) - 1)!! Binomial[2 p - 1 - j,
2 (p - j)] Binomial[j,
q] Pochhammer[-((alpha + m1)/2) - (j - q) + 1, j - q], {j, q,
p}], {q, 0, M}, {p, q, M}];
(*The coefficients Cp are polynomials in nu of order n-2 j for l=0 \
and 2 Floor(n-1)/2] - 2 j for l=1*)
(*Now we have: Cp[[1+l,1+n,1+Floor[(n-l)/2]]] \[Equal] ??*)

Cp = Table[
1 Sum[(-1)^(k + j - 1) Binomial[n,
q] ((( 2^(q) (q - k - 1/2)!) Pochhammer[(alpha + m1)/
2 - (n - q) + 1, n - q] )/(
Sqrt[Pi] Pochhammer[k + -j + (alpha + m1)/2,
2 j + l - k] )) Binomial[q, 2 (q - k)] Binomial[k - j - l,
j] , {q, j + l, n}, {k, Max[j + l, Ceiling[q/2]], q}], {l, 0,
1}, {n, 0, M}, {j, 0, Floor[(n - l)/2]}];

(*The coefficients AA[[1+m,1+j,1+k1,1+k2]] are ?.*)

MM = 10; Alpha =.; xi =.; nu =.;
d1[mm_, jj_] :=
Which[jj == 0,
2 (2 + Floor[1/4 (-3 + mm)] + Floor[1/4 (-2 + mm)] + Floor[mm/4]),
jj == 1,
Max[ 2 (1 + Floor[1/4 (-2 + mm)] + Floor[1/4 (-1 + mm)] +
Floor[mm/4]), 0]];
d2[mm_, jj_] :=
Which[jj == 0, 1 + Floor[1/2 (-1 + mm)] + 3 Floor[mm/2], jj == 1,
Max[1 + 3 Floor[1/2 (-1 + mm)] + Floor[mm/2], 0]];
AA = Table[
0, {mm, 0, MM}, {jj, 0, 1}, {kk1, 0, d1[mm, jj]}, {kk2, 0,
d2[mm, jj]}];
AA[[1, 1, 1, 1]] = 1;
For[mm = 0, mm <= MM - 1, mm++,
For[jj = 0, jj <= 1, jj++,
For[kk1 = 0, kk1 <= d1[mm + 1, jj], kk1++,
For[kk2 = 0, kk2 <= d2[mm + 1, jj], kk2++,
AA[[2 + mm, 1 + jj, 1 + kk1, 1 + kk2]] =
((Alpha + jj) - 2 mm + If[jj == 0, -nu, nu - 1]) If[
0 <= kk1 <= d1[mm, jj] + 0 && 1 <= kk2 <= d2[mm, jj] + 1,
AA[[1 + mm, 1 + jj, 1 + kk1, 0 + kk2]], 0]      +
(kk2 + 1) If[
2 <= kk1 <= d1[mm, jj] + 2 && 0 <= kk2 <= d2[mm, jj] - 1,
AA[[1 + mm, 1 + jj, -1 + kk1, 2 + kk2]] , 0]    +
(kk2 - 1) If[
0 <= kk1 <= d1[mm, jj] + 0 &&
2 <= kk2 <= d2[mm, jj] + 1,
AA[[1 + mm, 1 + jj, +1 + kk1, 0 + kk2]] , 0]  +

If[jj == 1 && 0 <= kk1 <= d1[mm, jj - 1] &&
1 <= kk2 <= d2[mm, jj - 1] + 1,
AA[[1 + mm, 0 + jj, 1 + kk1, 0 + kk2]], 0]  +

If[jj == 0 && 2 <= kk1 <= d1[mm, jj + 1] + 2 &&
1 <= kk2 <= d2[mm, jj + 1] + 1, -AA[[1 + mm,
2 + jj, -1 + kk1, 0 + kk2]], 0]  +

If[jj == 0 && 0 <= kk1 <= d1[mm, jj + 1] &&
3 <= kk2 <= d2[mm, jj + 1] + 3, -AA[[1 + mm, 2 + jj,
1 + kk1, -2 + kk2]], 0] ;
];
];
];
];

ll = (-I)^m2 2^(1/2 (-alpha + m1)) Sqrt[\[Pi]]  Table[
Binomial[(m1 - alpha - 1)/2, p] Pochhammer[
alpha + p + 1, (m1 - alpha - 1)/2 - p] CC[[1 + q,
p - q + 1]] x^(2 p + alpha)/2^p (
Sum[(-1)^(l - 1) Cp[[1 + l, 1 + q, 1 + j]] 1/2^(
2 j + l) (AA[[1 + m2, 1 + jj, 1 + 2 kk1,
1 + (m2 + 2 kk2 - 2 kk1)]] /. {Alpha :>
2 j + l - ((alpha + m1)/2) - 2 p,
nu :> (alpha + m1)/2 - l})
x^(2 kk1) y^(m2 + 2 kk2 - 2 kk1) (x^2 +
y^2)^((2 j + l - ((alpha + m1)/2) - 2 p + jj)/2 - m2)
BesselJ[(alpha + m1)/2 - l - jj, Sqrt[(x^2 + y^2)]], {l, 0,
1}, {j, 0, Floor[(q - l)/2]}, {jj, 0, 1}, {kk1, 0,
Floor[d1[m2, jj]/2]}, {kk2, 0, Floor[(m2 - jj)/2]}]), {q, 0,
M}, {p, q, M}];
N[Take[Accumulate[Flatten[ll]], -5], 15] // MatrixForm