# Skellam distribution: Deep connection between Poisson distributions and Bessel function?

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?

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. 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$).