Consider the generating function (here $z\in\mathbb{C}$, $|z|$ sufficiently small)
\begin{align*}
F(z)=\sum_{n=0}^{\infty}F_n z^n&:=\sum_{n=0}^{\infty}z^n\sum_{k=0}^{n}\frac{(k-n)^k}{k!}e^{n-k}
\\\color{gray}{[\text{replacing }k\text{ with }n-k]}\quad
&=\sum_{n=0}^{\infty}z^n\sum_{k=0}^{n}e^k\frac{(-k)^{n-k}}{(n-k)!}
\\\color{gray}{[\text{switching summations}]}\quad
&=\sum_{k=0}^{\infty}\sum_{n=k}^{\infty}z^n e^k\frac{(-k)^{n-k}}{(n-k)!}
\\\color{gray}{[\text{replacing $n$ with $n+k$}]}\quad
&=\sum_{k=0}^{\infty}(ez)^k\sum_{n=0}^{\infty}\frac{(-kz)^n}{n!}
\\\color{gray}{[\text{evaluating known sums}]}\quad
&=\sum_{k=0}^{\infty}(ze^{1-z})^k=\frac{1}{1-ze^{1-z}}.
\end{align*}
$F(z)$ has a double pole at $z=1$, with Laurent expansion $F(z)=2(z-1)^{-2}+(4/3)(z-1)^{-1}+\ldots$, and a sequence of simple poles (with the smallest in absolute value at $z\approx 3.0888\pm7.4615\mathrm{i}$). Thus, $$F(z)-\frac{2}{(1-z)^2}+\frac{4/3}{1-z}=\sum_{n=0}^{\infty}(F_n-2n-2/3)z^n$$ is regular in $|z|<r$ for some $r>1$ (we may take $r=8$ from the numerical value above).
In particular, the last series converges at $z=1$, which implies $\color{blue}{\lim\limits_{n\to\infty}(F_n-2n-2/3)=0}$.
In response to the comment by Gottfried Helms, here is a sketch of the analysis of $$y_k(x)=\sum_{n=0}^{\infty}\left\langle\begin{matrix}n\\k\end{matrix}\right\rangle\frac{x^n}{n!}\qquad(x\neq 0)$$ involving Eulerian numbers (denoted by $A(n,k)$ there). We consider the known $$Y(x,z):=\sum_{k=0}^{\infty}y_k(x)z^k=\frac{1-z}{e^{x(z-1)}-z}$$ as a function of $z$ (that is, with $x$ fixed); the denominator vanishes when $$(-xz)e^{-xz}=-xe^{-x}\iff z=z_m(x):=-W_m(-xe^{-x})/x\qquad(m\in\mathbb{Z})$$ where $W_m$ denotes the $m$-th branch of the Lambert W-function. Thus, $Y(x,z)$ has simple poles at these points, excluding $z=1$ if $x\neq 1$.
Now, an alternative representation of $y_k(x)$ comes from the partial fraction expansion of $Y(x,z)$, obtained using an approach going back to Cauchy, applicable to a meromorphic function $f(z)$ such that there is a positive integer $p$ and a sequence $\{C_n\}$ of simple contours, containing arbitrarily large circles inside, and $$\lim_{n\to\infty}\sup_{z\in C_n}|z^{-p}f(z)|=0.$$
In the case when $f(z)$ has only simple poles $z=z_m\neq 0$ with residues $a_m$, this approach yields $$f(z)=\sum_{k=0}^{p-1}\frac{f^{(k)}(0)}{k!}z^k+\sum_m\frac{a_m(z/z_m)^p}{z-z_m}.$$ For $f(z)=Y(x,z)$, the premise holds with $p=1$; computing $a_m$, we obtain $$y_k(x)=\sum_{m\in\mathbb{Z}}\frac{z_m(x)-1}{xz_m(x)-1}\big(z_m(x)\big)^{-k-1}\qquad(k>0)$$ which can be used to check the "empirical" results in the article being discussed.
No, the limit cannot exist. Say $\lim_{x\to +\infty} f(x) = L$. Name the sum $g: [0,\infty) \to \mathbb{R}$:
$$ g(x) = \sum_{k=0}^\infty \left[\lambda f\!\left(\frac{x}{2^k}\right) - f\!\left(\frac{x}{3 \cdot 2^k}\right)\right] $$
and suppose $\lim_{x \to +\infty} g(x) = M$.
Then for every positive $x$,
$$ g(2x) = \lambda f(2x) - f\!\left(\frac{2x}{3}\right) + g(x) $$
And then
$$ 0 = \lim_{x \to +\infty} \left[g(2x) - g(x) - \lambda f(2x) + f\!\left(\frac{2x}{3}\right)\right] = M - M + (\lambda - 1) L $$
So either $L=0$ or $\lambda=1$. It is not possible that $L \neq 0$ and $\lambda \neq 1$ and the limit of $g$ exists.
Best Answer
Your proof is nice, we have indeed
$$\left(1+\frac{1}{n}\right)^{n^2}\cdot\left(1+\frac{1}{n+1}\right)^{-(n+1)^2}=\\=\left(1+\frac{1}{n}\right)^{n^2}\cdot\left(1+\frac{1}{n+1}\right)^{-n^2}\cdot\left(1+\frac{1}{n+1}\right)^{-2n}\cdot\left(1+\frac{1}{n+1}\right)^{-1}=$$
$$=\left(\frac{(n+1)^2}{n(n+2)}\right)^{n^2}\cdot\left(1+\frac{1}{n+1}\right)^{-2n}\cdot\left(1+\frac{1}{n+1}\right)^{-1}=$$
$$=\left[\left(1+\frac{1}{n^2+2n}\right)^{n^2+2n}\right]^{\frac{n^2}{n^2+2n}}\cdot\left(1+\frac{1}{n+1}\right)^{-2n}\cdot\left(1+\frac{1}{n+1}\right)^{-1}\to e^1 \cdot e^{-2}\cdot 1 =\frac1e$$