I am interested in analytical formulas for integrals of products of spherical Bessel functions times a power:
$$I_{n,\ell}(r,R) \equiv \int_0^{\infty} dq \ q^{n} j_{\ell} \left( q r \right) j_{\ell}\left( q R \right) $$
where $n \in \mathbb{Z}$, $\ell \in \mathbb{N}$. This integral does not converge for every $n,l$. In that case I am interested in expressing the result as a distribution in terms of Dirac delta functions. I am interested in particular in the cases $n=-2,-1,1,3$.
Some results are already known, in particular for $n\ge 0$. For example:
$$\int \!\!dq j_{\ell} \left( q r \right) j_{\ell}\left( q R \right) =\frac{\pi}{2 \left( 1 +2 \ell\right)} \left( R^{-1-\ell} r^{\ell} \Theta_H \left( R- r \right) + r^{-1- \ell} R^{\ell} \Theta_H \left( r - R \right) \right)$$
where $\Theta_H$ denotes the Heaveside function. And for $n\ge 1$,
$$\int dq \ q^{2n} j_{\ell} \left( q r \right) j_{\ell} \left( q R \right) = \frac{\pi}{2 r^2} D^{(n-1)}_\ell(R)\delta_D \left( R - r \right)$$ where I have defined the differential operator acting on the Dirac delta by
$$D_\ell(R)\equiv-\frac{\partial^2}{\partial R^2} -\frac{2} {R}\frac{\partial}{\partial R} +\frac{\ell\left( \ell +1\right) }{R^2}$$
and $D^{(n-1)}_\ell(R)$ means that it is applied $n-1$ times. You can get these identities by using the differential equation satisfied by spherical Bessel functions.