The probability mass function for the [Skellam distribution](http://en.wikipedia.org/wiki/Skellam_distribution) for a count difference $k=n\_1-n\_2$ from two [Poisson-distributed](http://en.wikipedia.org/wiki/Poisson_distribution) 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](http://en.wikipedia.org/wiki/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?



Skellam distribution: Deep connection between Poisson distributions and Bessel function?
A mathematical reason is as follows.


On the one hand, the Laurent series for the modified Bessel functions of the first kind $I\_k$ can be deduced from the Laurent series for the Bessel functions of the first kind $J\_k$ given [here](http://en.wikipedia.org/wiki/Bessel_function). It reads 
$$
\sum\_{k\in\mathbb{Z}}I\_k(x)t^k=\mathrm{e}^{(x/2)(t+1/t)}.
$$
On the other hand, the characteristic function of a Poisson random variable $Y$ with mean $\mu$ is $E(z^{Y})=\mathrm{e}^{-\mu(1-z)}$, at least for every complex number $z$ of modulus $1$. Hence the characterization, which you recalled in your post, of Skellam distribution as the distribution of $X=Y\_1-Y\_2$ for independent Poisson random variables $Y\_1$ and $Y\_2$ with means $\mu\_1$ and $\mu\_2$ shows that 
$$
E(z^X)=\mathrm{e}^{-(\mu\_1+\mu\_2)}\mathrm{e}^{\mu\_1 z+\mu\_2/z}.
$$
Solving $(x/2)(t+1/t)=\mu\_1 z+\mu\_2/z$ for $x$ fixed and $t$ depending on $z$ yields $$
x=2\sqrt{\mu\_1\mu\_2},\qquad
t=z\sqrt{\mu\_1/\mu\_2}.
$$
The value of $\mathbb{P}(X=k)$ for every integer $k$ follows, which involves $I\_k(2\sqrt{\mu\_1\mu\_2})$. 


Finally, one should not worry too much about the appearance of $I\_k$ in this answer versus $I\_{|k|}$ in the OP's post because $J\_{-k}(-x)=(-1)^kJ\_k(x)$ for every integer $k$, hence $I\_{-k}(x)=I\_k(x)$ (a relation which is also a consequence of the invariance by $t\to1/t$ of the Laurent series for the functions $I\_k$).


