The probability mass function for the Skellam distribution for a count difference $k=n_1-n_2$ from two Poisson-distributed variables with means $\mu_1$ and $\mu_2$ is given by:

$$ f(k;\mu_1,\mu_2)= e^{-(\mu_1+\mu_2)} \left({\mu_1\over\mu_2}\right)^{k/2}I_{|k|}(2\sqrt{\mu_1\mu_2}) $$

where $I_k(z)$ is the modified Bessel function of the first kind.

**My question**: Is it just for convenience to get to grips with the resulting infinite summation terms that the Bessel function appears in this formula or is there a deeper mathematical reason connecting Poisson distributions and Bessel functions or even the Poisson distribution with Bessel differential equations? Is there perhaps even some physical interpretation or intution?