Roots of reverse Bessel polynomials define asymptotically a lattice

Question

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

Answer 1

Regarding asymptotics of real zero $\alpha_n$ ($2/r_{2n-1}$ in OP's notation) of generalized Bessel polynomials (with an extra parameter $a$ specialized in the current case to $a=2$), see M.G. de Bruin, E.B. Saff’ and R.S. Varga, On the zeros of generalized Bessel polynomials, Proceedings A, 84 (l), 20, (1981). In theorem $7.3$ they state

Obviously, $2/\alpha_n$ will be the real zero of the inverse polynomial. Here $\hat{r}=\lambda$ is Laplace limit in OP's notation, also changing $n\to 2n-1$ (since the numeration used in the paper is different), and noting $K(\hat{r};2)=\hat{r}$, we recover OP's formula.

Roots of Bessel polynomials asymptotically distribute along a certain curve in the complex plane, A. J. Carpenter, Asymptotics for the zeros of the generalized Bessel polynomials, Numer. Math. 62:465482 (1992).

To account for the zeros of the inverse polynomial one takes the conformal map $2/z$.