The free-particle propagator in 2D is given by $$ \rho_0(\mathbf{r},\mathbf{r'}, \beta) = \langle \mathbf{r'} \rvert e^{-\beta \hat{H}} \lvert \mathbf{r} \rangle = \frac{1}{4\pi \beta}e^{-\frac{(\mathbf{r}-\mathbf{r'})^2}{4\beta}}$$ where $\hat{H} = -\nabla^2_\mathbf{r}$ is the free-particle Hamiltonian in 2D, and where we assumed $\frac{\hbar^2}{2m} = 1$ and $\beta = \frac{it}{\hbar}$. This expression can easily be derived in Cartesian coordinates, cf my Phys.SE question here.
$$ $$
However, I am facing a two-particle scattering problem, first without interactions. In the center-of-mass frame, the Hamiltonian is the same as the one above with a reduced mass and I am trying to use polar coordinates (since the system has rotational symmetries).
Separating variables and taking only positive-energy solutions with $E=k^2\gt 0$, the wave-function can be written as $\psi(r, \theta) = \Theta_l(\theta)R_{kl}(r)$. The Schrödinger equation yields $\Theta_l(\theta) = \frac{1}{\sqrt{2\pi}}e^{il\theta}$ and $R_{kl}$ is a solution to
$$ \frac{\mathrm{d}^2R}{\mathrm{d}r^2} + \frac{1}{r}\frac{\mathrm{d}R}{\mathrm{d}r} + (k^2-l^2)R=0.$$
A general solution to this equation is given by
$$ R_{kl}(r) = \mathcal{C}_k [\cos\delta_k J_l(kr) + \sin\delta_kY_l(kr)]$$ where $J_l$ and $Y_l$ are Bessel functions of the first and second kind, $\mathcal{C}_k$ is a normalization constant and $\delta_k$ is the phase shift to be determined.
Since there are no interactions between the two particles, the solutions in $Y_l$ (which diverge at $r=0$) should vanish, and therefore $\delta_k = 0$. Furthermore, since $R_{kl}$ is in units of $(\text{distance})^{-1/2}$, we take $\mathcal{C}_k = \sqrt{k}$.
We can now derive the propagator with
$$ \begin{align}\rho_0(\mathbf{r},\mathbf{r'}, \beta) & = \sum_l \sum_k \psi_{k,l}^*(\mathbf{r'}) e^{-\beta E} \psi_{k,l}(\mathbf{r}) \\ & = \frac{1}{2\pi} \sum_l e^{il(\theta-\theta')}\int_0^\infty k J_l(kr')e^{-\beta k^2} J_l(kr) \mathrm{d}k \\ \rho_0(\mathbf{r},\mathbf{r'}, \beta) & = \frac{1}{4\pi\beta} e^{-\frac{r^2 + r'^2}{4\beta}} \sum_l e^{il(\theta-\theta')} I_l\left(\frac{rr'}{2\beta}\right) \end{align}$$
where $I_l$ is the modified Bessel function of the first kind. This last integral is given in DLMF. This very much resembles the first equation in the post, but Mathematica says the two are not equal. How can I reconcile the two results?
