(Title changed after comment from Nemo).

Set $$P_n=\sum_{k=0}^n\frac{(n+k)!x^{n-k}}{k!(n-k)!}\ .$$ (The polynomial $\theta_n(x)=2^{-n}P_n(2x)$ is the so-called reverse Bessel polynomial, see comment from Nemo.)

Let $\mathcal R_n$ denote the set of all roots of $P_n$. The union $\bigcup_{n=1}^\infty \mathcal R_n$ intersects a small window centered at a very large negative real number almost in a sort of limit-lattice $\Lambda$ . (The sets $\mathcal R_n$ are close to parabolas opening to the right.) The real generator of $\Lambda$ (defined as $\lim_{n\rightarrow\infty}r_{2n-1}-r_{2n+1}$ where $r_{2n-1}=\mathcal R_{2n-1}\cap \mathbb R$ is the real root of $P_{2n-1}$) seems to be given by $4\lambda$ where $$\lambda=0.662743419349\ldots$$ is the real solution of $$\lambda e^{\sqrt{1+\lambda^2}}=1+\sqrt{1+\lambda^2}\ .$$ (The constant $\lambda$, sometimes called "Lagrange's constant", is apparently related to celestial mechanics.)

One seems to have $$\lim_{n\rightarrow \infty}\frac{r_{2n-1}}{4n-1}=-\lambda\ .$$ Experimentally, there exists seemingly a sequence \begin{align*} a_2&=0.05490008226\ldots,\\ a_4&=-0.01648163\ldots,\\ a_6&=0.021761\ldots,\\ &\ \ \vdots \end{align*} such that we have asymptotic expansions $$-\frac{r_{2n-1}}{(4n-1)\lambda}=1+\sum_{n=1}^N\frac{a_{2n}}{(4n-1)^{2n}}+O(n^{-2N-2})$$ involving the unique real root $r_{2n-1}=\mathcal R_{2n-1}\cap\mathbb R$ of $P_{2n-1}$.

I noticed nothing special for "the" other generator of $\Lambda$.

Are there any explanations for this curious behaviour?

(Added after comment by Nemo: Since Bessel polynomials are orthogonal polynomials for a suitable measure on the complex unit-circle and are useful for filters, an explanation involves probably properties of Bessel polynomials.)