EDITADO 2. Este es otro enfoque: El punto de partida es la siguiente identidad
$$ \int_{0}^{\infty} \frac{(1 - e^{-x})^{n}}{x^{m}} \, dx = \frac{(-1)^{m}}{(m-1)!} \sum_{k=1}^{n} \binom{n}{k} (-1)^{k} k^{m-1} \log k. \tag{1} $$
que probé aquí . Enchufando $m = n$ a $\text{(1)}$ tenemos
$$ A_{n} := \frac{(-1)^{n}n^{2}}{n!} \sum_{k=2}^{n} \binom{n}{k}(-1)^{k} k^{n-1} \log k = \int_{0}^{\infty} n \left( \frac{1 - e^{-x}}{x} \right)^{n} \, dx. $$
Ahora dejemos que $f(x) : [0, \infty) \to (0, 1]$ se define por $f(x) = (1 - e^{-x})/x$ . Esta función es monótona y decreciente, y denotamos su inversa por $g = f^{-1}$ . Entonces se deduce de la sustitución $x = g(y)$ que
$$ A_{n} = \int_{0}^{1} n y^{n} \left| g'(y) \right| \, dy. $$
Al observar que $y \mapsto ny^{n}$ se comporta como una aproximación a la identidad de la masa unitaria $\delta_{1}$ se deduce que
$$ \lim_{n\to\infty} A_{n} = \left| g'(1) \right| = \frac{1}{\left| f'(0) \right|} = 2. $$
EDITADO. Por fin he dado con una solución completa:
Para no tener en cuenta la cancelación entre grandes cantidades, encontramos una representación más manejable. Sea
$$ A_{n} = \frac{(-1)^{n}n^{2}}{n!} \sum_{k=2}^{n} \binom{n}{k}(-1)^{k} k^{n-1} \log k $$
e introducir
$$ f_{n}(x) = \frac{n^{2} x^{n-2} \log x}{(x-1)\cdots(x-n)}. $$
Entonces tenemos
$$ \operatorname{Res}\limits_{z=k} \, f_{n}(z) = \frac{(-1)^{n-k} n^{2} k^{n-1} \log k}{k!(n-k)!} $$
y, por tanto, para cualquier contorno cerrado simple $C \subset \Bbb{C} \setminus (-\infty, 0]$ bobinado $\{1, \cdots ,n\}$ en sentido contrario a las agujas del reloj se deduce que
$$ \int_{C} f_{n}(z) \, dz = 2\pi i A_{n}. $$
Pero como $f_{n}(z) = O(|z|^{-2}\log|z|)$ como $|z| \to \infty$ , podemos modificar el contorno de manera que $C$ vientos $(-\infty, 0]$ y un simple cálculo muestra que
$$ A_{n} = \int_{0}^{\infty} \frac{n^{2} x^{n-2}}{(x+1)\cdots(x+n)} \, dx. \tag{1} $$
Utilizando la sustitución $x \mapsto n^{2}/x$ se deduce de $\text{(1)}$ que
$$ A_{n} = \int_{0}^{\infty} \prod_{k=1}^{n} \left( 1 + \frac{kx}{n^{2}} \right)^{-1} \, dx. \tag{2} $$
Ahora bien, observe que cuando $n \geq 2$ expandiendo el denominador y descartando el término de orden superior, tenemos
\begin{align*} \prod_{k=1}^{n} \left( 1 + \frac{kx}{n^{2}} \right) &\geq 1 + \sum_{1\leq k \leq n} \frac{kx}{n^{2}} + \sum_{1 \leq l < k \leq n} \frac{kl x^{2}}{n^{4}} \\ &= 1 + \frac{n+1}{2n} x + \frac{(3n+2)(n^{2}-1)}{24n^{3}} x^{2} \\ &\geq 1 + \frac{x}{2} + \frac{x^{2}}{8}. \end{align*}
Esto significa que las integradas de $\text{(2)}$ están limitadas por una función integrable $1/(x^{2}/8 + x/2 + 1)$ y así podemos aplicar el teorema de convergencia dominada, una vez que probamos que el integrando converge puntualmente. Pero es inmediato que para cualquier $x$ ,
$$ \log \prod_{k=1}^{n} \left( 1 + \frac{kx}{n^{2}} \right) = \sum_{k=1}^{n} \log \left( 1 + \frac{kx}{n^{2}} \right) = \frac{x}{2} + \mathcal{O}\left(\frac{1}{n}\right) \longrightarrow \frac{x}{2} $$
y por lo tanto
$$ \lim_{n\to\infty} A_{n} = \int_{0}^{\infty} e^{-x/2} \, dx = 2. $$