Evaluating the sum $f(x):=\sum_{n=1}^\infty \frac{1}{n! n^n}(-x^2)^n$ and estimating bounds For real variable $x$, the function
\begin{equation}
f(x):=\sum_{n=1}^\infty \frac{1}{n! n^n}(-x^2)^n
\end{equation}
clearly has infinite radius of convergence and defines a $C^\infty$ function on $\mathbb{R}$.
However, I wonder if this $f(x)$ is a bounded function on $\mathbb{R}$ as well. Also, is it possible to estimate the bound if it is indeed bounded?
I tried to apply summation by parts or the error bounds for an alternating series. But all I have found out is that the partial sums $\sum_{n=1}^N \frac{1}{n!}(-x^2)^n$ are extremely slowly convergent to the Gaussian function..
Could anyone please help me with the above function $f(x)$?
 A: Due to the bounty I'll give a second answer.
Let $r = \sqrt{-8\log \tau}$.
The proposed asymptotic form of
$$
f(x)= - \int_0^1 \mathrm d\tau \,x  r \, J_1(x \tau r),\tag{1}
$$ for large $x$ from my first answer,
$$
f_\infty(x)= \sqrt 2 \cos(2 e^{-1/2} x) + c + o(1),\tag{2}
$$
(note that the $O(x^{-1/2})$ term was too optimistic) can be derived and refined in the following way:
Obviously $f_\infty(x)$ asymptotically ($\simeq$) fulfills the ODE (see also the comment of @TheSimpliFire)
$$
f_\infty(x)+\frac e 4 f_\infty''(x) \simeq c,\tag{3}
$$
therefore we check if $f(x)$ also solves (3) for large $x$.
We get
\begin{align}
F(x) &= f(x)+\frac e 4 f''(x) \tag{4a}\\
&=\int_0^1 \mathrm d\tau\,
\frac r 4\big[x(e r^2t^2-4) J_1(x\tau r)-e r t J_0(x\tau r)\big]\tag{4b}\\
&\simeq\int_0^1 \mathrm d\tau\,
\frac {xr} 4(e r^2t^2-4) J_1(x\tau r).\tag{4c}
\end{align}
Note that for large $x$ the integrand oscillates and decays rapidly with increasing $x$, such that the relevant contributions come from the lower bound, where $\tau\ll 1$.
If we now substitute $\tau\to y/x$ and use $-\log \tau \to \log x-\log y \simeq \log x$, we get
\begin{align}
F(x)&\simeq-\sqrt{8\log x}\int_0^\infty \mathrm dy\,J_1\big(\sqrt{8\log x}\,y\big) =-1,\tag{5}
\end{align}
such that asymptotically $c=-1$. Note that the convergence is extremely slow.
In conclusion, the asymptotics becomes
$$
f_\infty(x)= \sqrt 2 \cos(2 e^{-1/2} x) - 1 + o(1),\tag{6}
$$
the constant $c=-1$ can be seen as the missing contribution from the $n=0$ term in the original sum (if we define $0^0=1$).
The next terms in the expansions seem to be
$$
f_\infty(x) \simeq \sqrt 2 \cos(2 e^{-1/2} x) - 1 + d x^{-\epsilon} + e x^{-1}\sin(2 e^{-1/2} x),\tag{7}
$$
with $d\approx -0.18$, $\epsilon\approx 0.17$, $e\approx -0.097$.
A: There is some evil spell on this problem because it somehow makes many very smart and reputable people to write total nonsense as an answer and then delete it (and not so smart people like myself fall into the same trap: I already posted a nonsensical computation twice!).
Let me try to lift this spell now. Everybody agrees that the key is the integral representation in terms of Bessel functions. I just want to use as little high-tech as possible. My starting point is the Sophomore's dream whose underlying identity is
$$
\frac 1{(n-1)!}\int_0^1 a(x)^{n-1}\,dx=\frac 1{n^n}
$$
with $a(x)=-x\log x$.
This leads to the representation (with $t\ge 0$ having the meaning of the original $x^2$)
$$
\sum_{n\ge 1}\frac{(-1)^n}{n^n n!}t^n=t\int_0^1\left[\sum_{n\ge 1}\frac{(-1)^n}{(n-1)!n!}(a(x)t)^{n-1}\right]\,dx
\\
=t\int_0^1 G(a(x)t)\,dx
$$
where
$$
G(t)=\sum_{n\ge 1}\frac{(-1)^n}{(n-1)!n!}t^{n-1}=H'(t)
$$
with
$$
H(t)=\sum_{n\ge 0}\frac{(-1)^n}{(n!)^2}t^n=J_0(2\sqrt t)
$$
where $J_0$ is the classical Bessel function.
Assuming that the estimates $J_0(t),J_0'(t)=O(t^{-1/2})$ as $t\to +\infty$ are known, we conclude that $G(t)=O(t^{-3/4})$ and $H(t)=O(t^{-1/4})$ as $t\to +\infty$. Also both are uniformly bounded.
Note now that $a(x)$ attains its maximum at $x_0=e^{-1}$ and $a'(x)=-1-\log x$ is a monotone function preserving sign on each of the intervals $[0,x_0)$ and $(x_0,1]$. Now if we put $x_{\pm}=x_0\pm t^{-1/4}$ for not too small $t$ (the boundedness on compact intervals is trivial), then $t\int_{[x_-,x_+]}G(a(x)t)\,dx$ is bounded independently of $t$ due to the upper bound for $|G|$ and the fact that $a(x)$ is uniformly separated from $0$ on that interval.
On the remaining set, we can write $tG(a(x)t)=\frac 1{a'(t)}\frac d{dt}H(a(x)t)$. Now we can use a simple lemma from undergraduate calculus:
If $I$ is an interval, $U$ is a smooth bounded function on $I$ and $v$ is a monotone function on $I$, then
$$
\left|\int_I U'v\right|\le 4\|U\|_\infty \|v\|_\infty
$$
Applying it to $I=[0,\frac{x_0}2]$ and $I=[\frac{x_0}2,x_-$] with $U(x)=H(a(x)t)$ and $v(x)=\frac 1{a'(x)}$, we get the bounds ${\rm const}\cdot{\rm const}$ and $Ct^{-1/4}\times Ct^{1/4}$ respectively. The interval $[x_+,1]$ is treated in the same way.
It remains to give an elementary proof of the bounds for $G=H'$ and $H$ for those poor mortals like myself who can never remember fancy names and properties of special functions.
Let's look at $H$. It satisfies the differential equation $(tH')'+H=0$ (straight from the power series expansion). We'll now prove the following simple lemma (crude JWKB for those who like fancy names):
Let $A>0$ be a reasonably smooth and not too fast changing function on $[0,+\infty)$ (the exact conditions will be derived below) and let $u$ solve the differential equation $(Au')'+u=Au''+A'u'+u=0$. Then we can form the energy
$$
E=A(u')^2+u^2
$$
and the shmark (I'm sure there is a well-established name for it, but I am just an ignoramus in such matters)
$$
S=A(u')^2-u^2\,.
$$
We have
$$
E'=2Au'u''+A'(u')^2+2uu'=2u'(Au''+A'u'+u)-A'(u'^2)=-A'(u')^2
\\
=-\frac{A'}{2A}[1+\tfrac SE]E
$$
In particular, we see that $\frac{|E'|}{E}\le \frac{|A'|}{A}$.
Also if $W=Au'u$, then
$$
W'=(Au')'u+A(u')^2=-u^2+A(u')^2=S
$$
and, by Cauchy Schwarz, $|W|\le\sqrt A E$.
Now, from the differential equation for $E$, we have
$$
\log E(t)=\log E(0)-\int_0^t \frac{A'}{2A}(1+\frac SE)\,.
$$
Integrating $\frac{A'}{2A}$, we get $C-\frac 12\log A(t)$. Now,
$$
\int_0^t \frac{A'}{2A}\frac SE=\int_0^t W'\frac{A'}{2AE}=
W\frac{A'}{2AE}|_0^t-\int_0^t W\frac{A'}{2AE}[\frac {A''}{A'}-\frac{A'}A-\frac{E'}E]\,.
$$
We want to show that it is just bounded. For that we can use $|W|\le \sqrt AE$ and $\frac {|E'|}{E}\le \frac{|A'|}A$. Thus the set of sufficient conditions is
$$
\frac{|A'|}{\sqrt{A}}=O(1), \int^\infty\left[\frac{|A''|}{\sqrt{A}}+\frac{|A'|^2}{A^{3/2}}\right]<+\infty
$$
(all near $+\infty$, what happens near $0$ is of no importance).
Then we conclude that $E=O(A^{-1/2})$ and, thereby $u=O(A^{-1/4})$, $u'=O(A^{-3/4})$ at $+\infty$. It is easy to check that $A(t)=t$ satisfies the sufficient assumptions. I talked to Alexandre Eremenko and he said that he had always derived any JWKB related estimate by using iterations and analytic continuation in his classes, so I decided to post this elementary energy derivation that requires no analyticity and runs directly on the physical idea that the kinetic and the potential energies are balanced on average.
A: As already stated by @Noam and @Alexandre in the comments above,
$$
f(x)= - x \int_0^1 \mathrm d\tau \sqrt{-8\log \tau} \, J_1\big(x \tau \sqrt{-8\log \tau}\big),\tag{1}
$$
with Bessel function $J_1$. Here, I did one further substitution in order to simplify the expression. The function $f(x)$ oscillates at large $x$, and asymptotically seems to become
$$
f_\infty(x)= \sqrt 2 \cos(2 e^{-1/2} x) + c + O(x^{-1/2}),\tag{2}
$$
with $c\approx -1.055$. The oscillation frequency stems from the oscillations of  $J_1(x \tau \sqrt{-8\log \tau})$ in the integrand at stationary maximal argument (see comment by @Noam), which is at $\tau=e^{-1/2}=0.60653\ldots$. The prefactor of $x^{-1/2}$ is around $-0.285$.
