Solving an integral involving a Bessel function, Laguerre function and Gaussian We want to calculate the expectation value $\langle q^2\rangle$ in polar coordinates which gives us the following integral, for integer values of $p$:
\begin{equation}\int_0^\infty dq~q^3 \left(\int_0^\infty dr~ J_{0}(q r) e^{-r^2/2} L_{p}^{1}(r^2) r^2\right)^2.\end{equation}
With Bessel function of the first kind $J_n(x)$ and generalized Laguerre polynomial $L_{p}^{l}(x)$.
We can't solve it because we don't find a nice form of the Fourier-Bessel transform of our waves: (with $l\geq1$)
\begin{equation}\psi_p^l(q)=\int_0^\infty dr~ J_{l-1}(q r) e^{-r^2/2} L_{p}^{l}(r^2)r^{l+1}.\end{equation}
In the literature (http://hdl.handle.net/1794/3778, Eq. 3.26) we find a form which is similar but not quite the same.
\begin{equation}\int_0^\infty dr~ J_{l}(q r) e^{-r^2/2} L_{p}^{l}(r^2)r^{l+1}=(-1)^p q^l L_p^l(q^2) e^{-q^2/2}.\end{equation}
Can anyone help us to solve the first or second integral?
 A: This is my solution. Might still contain mistakes. Hope this helps others too!
The bessel function
Using the following properties
\begin{equation}
    \int_0^\infty dq q J_v(q \rho_1) J_v(q \rho_2)= \frac{\delta(r_1-r_2)}{r_1}
\end{equation}
\begin{equation}
    J'_v(z)=- J_{v+1}(z)+\frac{v}{z}J_{v}(z)
\end{equation}
\begin{equation}
    J_0'(z)=-J_1(z)
\end{equation}
\begin{equation}
    J_{-n}=(-1)^n J_{n}.
\end{equation}
Therefore;
\begin{equation}
    q^3 J'_{-1}(q \rho_1) J'_{-1}(q \rho_2)= q^3 J_0(q \rho_1) J_0(q \rho_2)+ q \frac{1}{\rho_1 \rho_2} J_{-1}(q\rho_1)J_{-1}(q\rho_2)+q^2\left( \frac{1}{\rho_2}J_0(q \rho_1) J'_0(q \rho_2)+\frac{1}{\rho_1}J'_0(q \rho_1) J_0(q \rho_2) \right)
\end{equation}
This gives
\begin{equation}
    \int_0^\infty dq q^3 J_0(q \rho_1) J_0(q \rho_2)=\left(\partial_{\rho_1}\partial_{\rho_2}-\frac{1}{\rho_1 \rho_2}-\frac{\partial_{\rho_1}}{\rho_1}-\frac{\partial_{\rho_2}}{\rho_2}\right)\frac{\delta(\rho_1-\rho_2)}{\rho_1}
\end{equation}
Solving the integral
\begin{equation}
    I=\int_0^\infty d \rho f(\rho) \left(-\frac{\partial^2_{\rho}}{\rho}+\frac{3\partial_{\rho}}{\rho^2}-\frac{3}{\rho^3} \right)f(\rho).
\end{equation}
define $x\equiv \rho^2$.
\begin{equation}
    I=-\int_0^\infty dx f(x)\left[ 2\partial^2_x-2\frac{\partial_x}{x}+\frac{3}{2x^2}\right] f(x). 
\end{equation}
Using the identities
\begin{equation}
    \partial_x(x L^1_p(x))=(p+1) L_{p}^0
\end{equation}
\begin{equation}
    \partial_x(L_p^0)=-L_{p-1} (\text{for}~ p\geq1) ~(=0~ \text{otherwise})
\end{equation}
the derivatives to $f(x)=x L_p^1(x) e^{-x/2}$ are
\begin{equation}
\partial_x f(x)=(- \frac{x}{2} L_p^1(x)+(p+1)L_p^0(x) ) e^{x/2}
\end{equation}
\begin{equation}
\partial^2_x f(x)=(\frac{x}{4} L_p^1(x)-(p+1)L_p^0(x)-L^1_{p-1}(!)) e^{x/2}
\end{equation}
where the $(!)$ points out that this term is only there is $p\geq1$.
This gives the following integral
\begin{equation}
    I=-\int_0^\infty dx e^{-x} \left[L_p^1(x)L_p^1(x)(\frac{3}{2}+x+\frac{x^2}{2}) -x (p+1) L_p^0(x) L_p^1(x)- 2 x L_{p-1}^1(x) L_{p}^1(x)\right]
\end{equation}
Orthogonality conditions
Making us of the orthogonality conditions
\begin{equation}
    \int_0^\infty dx x^a e^{-x} L_n^a L_m^a =\frac{(n+a)!}{n!} \delta_{n,m}
\end{equation}
\begin{equation}
    \int_0^\infty dx x^{a+1} e^{-x} (L_n^a(x))^2 =\frac{(n+a)!}{n!} (2n+a+1)
\end{equation}
and
\begin{equation}
    L_p^1=\sum_{i=0}^p L_{i}^0(x)
\end{equation}
\begin{equation}
    L_n^a=L_n^{a+1}-L_{n-1}^{a+1}.
\end{equation}
then
\begin{equation}
    I=-(3p/2+(p+1)+(p+1)(2p+1+1)/2-(p+1)(p+1)-0)=-\frac{5p}{2}-1
\end{equation}
