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

Question

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?

Answer 1

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