I have asked similar questions on Math Stack Exchange, but not been able to receive many helpful responses. Therefore, I am posting this problem here, and any input would be extremely valuable.

I have the following integral:

$$\displaystyle \int_{|k| \geqslant 1} \left( \int_{|\alpha| \geqslant 1} \frac{J_{1}(\rho |\alpha|)J_{1}(\rho|k- \alpha|)}{|\alpha||k-\alpha|} \ \mathrm{d}\alpha \right)^2 \ \mathrm{d}k,$$

where both $\alpha$ and $k$ are vectors in $\mathbb{R}^2$, with $k \neq 0$, $\rho \gg 1$ does not depend on either $k$ nor $\alpha$, and $J_{\nu}$ denotes the Bessel function of the first kind. I'm having some trouble with the best way to approach this integral. If we focus on the inner integral first, then using the fact that for sufficiently large, positive $z$ we have $|J_{\nu}(z)| \leqslant C|z|^{-1/2},$ then the inner integral can be reduced to

$$\displaystyle \int_{|\alpha| \geqslant 1} |\alpha|^{-3/2}|k-\alpha|^{-3/2} \ \mathrm{d}\alpha.$$

However, as can be seen in this answer, this integral is $O(|k|^{-1})$, which, after squaring, is clearly not integrable over all $|k| \geqslant 1$ after switching to polar co-ordinates. We would need an estimate of at least $O(|k|^{-1 - \epsilon})$ for any $\epsilon > 0$ to guarantee convergence of the outer integral.

One idea might be to try to bring the outer integral inside (though one would need to justify interchanging the order of integration). Using the asymptotics for the Bessel functions gives a product of cosines, and then one can use polar co-ordinates (taking $r = |\alpha|$). This would cancel out the $|\alpha|$ in the denominator, but then the $|k-\alpha|$ terms get very messy, and this might make things worse. The Bessel functions appear to cause the most trouble. Does anyone have any ideas on how to proceed?

This integral arises when considering a problem on the distribution of lattice points inside a ball in $\mathbb{R}^d$ with radius $\rho$, which can be thought of as a generalisation of some aspects of the Gauss circle problem. For the full details, I recommend this paper (particularly pages 10-11 and 15-16). The paper considers $\sigma_p$ and proves some asymptotic bounds for $p = 1$ and $p = 2$ which depend on the dimension $d$ and the radius $\rho$. I'm attempting to generalise that work to $p = 4$, and with some extra work, any $p \in \mathbb{N}.$ The derivation of this integral begins here, and uses the estimate at the bottom of this post to arrive at the aforementioned integral.

My expectation is that the entire expression should be $O(\rho^{-2})$. (This is the information that I would like to extract from the integral.) This would be the best possible bound achievable, as it is known that the integral has the same power of $\rho$ as a lower bound.