If $X$ and $Y$ are two independent Poisson random variables with same parameter $\lambda$, $Z=X-Y$ is the Skellam distribution with parameters $\lambda$ and $\lambda$. This is distributed around zero with a width which scales with the value of $\lambda$.
I found numerically that the distribution of $\frac{X-Y}{\sqrt{+}}$ is independent of $\lambda$. Does someone know if this PDF has a known analytical form, or how I could try to derive it? Also interested in the normal approximation if its easier.
I will here be giving a normal approximation, which will hold for large $\lambda$. Using the central limit theorem we can conclude, that $$ \sqrt{\lambda}(\frac{1}{\lambda}\begin{pmatrix} X \\\ Y \end{pmatrix}-\begin{pmatrix} 1 \\\ 1 \end{pmatrix}) \xrightarrow{\sim} N_2(0, \begin{pmatrix} 1 & 0 \\\ 0 & 1 \end{pmatrix}) \quad \text{as $\lambda \rightarrow \infty$}. $$ If we let $g(x,y)=\frac{x-y}{\sqrt{x+y}}$, then the asymptotic distribution of $g(X/\lambda,Y/\lambda) = \frac{1}{\sqrt{\lambda}}\frac{X-Y}{\sqrt{X+Y}}$ can be found by the delta method. The delta method gives us that $$\frac{X-Y}{\sqrt{X+Y}} = \sqrt{\lambda}(g(X/\lambda , Y \lambda)-g(1,1)) \xrightarrow{\sim} N(0, || \nabla g(1,1) ||^2).$$ Indeed we can see that the asymptotic distribution does not depend on $\lambda$. Let us calculate the asymptotic variance. We note that $$\nabla g(x,y) = \begin{pmatrix} \frac{x+3y}{2(x+y)^{(3/2)}} & \frac{-3x-y} {2(x+y)^{(3/2)}})\end{pmatrix},$$ yielding that $$||\nabla g(1,1)||^2 = \big(\frac{1+3}{2(1+1)^{(3/2)}}\big)^2 + \big(\frac{-3-1}{2(1+1)^{(3/2)}}\big)^2 = 1,$$ We may thus for large values of lambda conclude that $\frac{X-Y}{\sqrt{X+Y}} \sim N(0,1)$ approximately.