Quiero añadir un comentario a Noam D. Elkies' hermosa respuesta. A partir de la representación integral para $f$, poniendo a $e^{-t}=s$ en la integral,
$$f(-x)=1-x\int_0^\infty e^{-xte^{-t}} e^{-t}dt = 1-x\int_0^1 s^{sx}ds\, ,$$
así que, para $x\to \infty$, $ f(-x)=o(1)$ es equivalente a $$\int_0^1 xu(s)^xds=1+o(1)\, ,$$ where $u\en C([0,1])$ is the function $u(s):=s^s$. As a matter of fact, since $0\le u(s)\le 1$ for all $s$ and $u(s)=1$ only for $s=0$ or $s=1$, it turns out that the limit only depends on $u'(0)$ and $u'(1)$.
Desde $u'(1)=1$, para cualquier $\lambda < 1 < \mu$ existe un $b < 1$ tal que
para todos los $s\in [b,1]$ no tiene
$$1+\mu(s-1) \le u(s)\le 1+\lambda(s-1)\, ,$$
así que
$$x\big(1+\mu(s-1)\big)^x \le xu(s)^x\le x\big(1+\lambda(s-1)\big)^x\, .$$
Del mismo modo, desde la $u'(0)=-\infty$, para cualquier $\nu > 0$ existe un $a > 0$ tal que
para todos $s\in [0,a]$ $$u(s)\le1-\nu s\, ,$$ así
$$xu(s)^x\le x\big( 1-\nu s\big) ^ x \, .$$
Por otra parte, ya que en cualquier intervalo de $[a,b]\subset\subset(0,1)$ la función de $u$ se apartó de $1$, está claro que $\int_a^b xu(s)^xds=o(1)$ por la convergencia uniforme a $0$.
La integración de más de $s\in [ 0,1]$, y recordando que $\lambda < 1 < \mu$ e $\nu > 0$ fueron arbitrarias, las desigualdades por encima claramente de dar
$$\int_0^1 xu(s)^xds=\int_0^a x u(s)^xds+\int_a^b xu(s)^xds+\int_b^1 xu(s)^xds=1+o(1) \, ,$$
para $x\to \infty$.