Closed-form or integral representation for $\sum_{n=0}^\infty\frac{e^{-\lambda}\lambda^n}{n!}\phi(z-n;\mu;\sigma)$. Let $Z=X+Y$, where $X\sim\operatorname{Poisson}(\lambda)$ and $Y\sim\mathcal N(\mu,\sigma^2)$ are independent.  The density of $Z$ can thus be expressed as an infinite component Gaussian mixture of the form
$$
f_Z(z)=\sum_{n=0}^\infty\frac{e^{-\lambda}\lambda^n}{n!}\phi(z-n;\mu;\sigma).
$$
Here, $\phi$ denotes the normal density with mean $\mu$ and variance $\sigma^2$.  My question is if $f_Z$ can be expressed in a "closed-form" that is compatible with numerical computation, e.g. in terms of known special functions or an integral.  I am also interested in any special cases where such a closed-form may be obtained, for example, when $\mu=0$.
A related question assigns equal weights to each Gaussian component which subsequently allows for an integral representation to be derived that is expressible in terms of error functions.
In terms of special functions the Jacobi theta function $\vartheta_3$ is the only I'm aware of that admits a series expansion in terms of Gaussian-like functions.
For the (perhaps) most basic of special cases: $\mu=0$, $\sigma=1$, and $\lambda=1$, I cannot obtain a closed-form and Mathematica returns an unevaluated series.
Any leads are also welcome.

Edit:
As pointed out by @K.defaoite we can write $f_Z$ as
$$
f_Z(z)=\phi(z;\mu,\sigma)e^{-\lambda}\sum_{n=0}^\infty\frac{x^n}{n!}e^{-n^2/(2\sigma^2)},
$$
where $x=\lambda\exp((z-\mu)/\sigma^2)$. This makes the series a little simpler to work with.
 A: An integral representation can be obtained from a result proposed by M.D. Schmidt in Square series generating function transformations (Prop. 5.2):
\begin{align}
 E_{\mathrm{Sq}}(q, r, z)&=\sum_{n=0}^\infty q^{n^2}r^n\frac{z^n}{n!}\\
  &=\int_{0}^{\infty} \frac{e^{-t^{2} / 2}}{\sqrt{2 \pi}}\left[e^{e^{t\sqrt{2 \log (q) }} r z}+e^{e^{-t\sqrt{2 \log (q)}} r z}\right]  \,dt
\end{align}
when $|q|<1,r,z\in \mathbb C$.
Here, we have
\begin{equation}
 \sum_{n=0}^\infty\frac{x^n}{n!}e^{-n^2/(2\sigma^2)}=E_{\mathrm{Sq}}(e^{-1/(2\sigma^2)}, 1, x)
\end{equation}
and thus,
\begin{equation}
 \sum_{n=0}^\infty\frac{x^n}{n!}e^{-n^2/(2\sigma^2)}=\sqrt{\frac{2}{\pi}}\int_0^\infty e^{-t^2/2+x\cos t/\sigma}\cos\left( x\sin\frac t\sigma \right)\,dt
\end{equation}
which seems to be numerically correct.
A: It occurred to me that this can be accomplished via characteristic functions. Since $X$ and $Y$ are independent
$$
\varphi_Z(t)=\mathsf Ee^{itZ}=\varphi_X(t)\varphi_Y(t).
$$
The characteristic functions of the Poisson and normal densities are known yielding
$$
\varphi_Z(t)=\exp(\lambda (e^{it}-1))\exp(i\mu t-\sigma^2t^2/2).
$$
To find an integral representation for the density $f_Z(z)$ we may apply the characteristic function inversion formula
$$
f_{Z}(z)=\frac{1}{2\pi}\int_{-\infty}^\infty e^{-itz}\varphi_{Z}(t)\,\mathrm dt.
$$
After a couple applications of Euler's identity, the integrand can be written as
$$
e^{-itz}\varphi_{Z}(t)=
\underbrace{e^{\lambda(\cos t-1)-\sigma^2t^2/2}\cos(\lambda\sin t+t(\mu-z))}_{\text{even}}
+i\underbrace{e^{\lambda (\cos t-1)-\sigma^2t^2/2}\sin(\lambda\sin t+t(\mu-z))}_{\text{odd}}
$$
The odd term does not contribute to $f_Z$. Likewise, the even term can be rewritten by multiplying the integrand by two and integrating over $[0,\infty)$ leaving us with
$$
f_{Z}(z)=\frac{1}{\pi}\int_0^\infty e^{\lambda(\cos t-1)-\sigma^2t^2/2}\cos(\lambda\sin t+t(\mu-z))\,\mathrm dt;
$$
thus,
$$
\sum_{n=0}^\infty\frac{e^{-\lambda}\lambda^n}{n!}\phi(z-n;\mu;\sigma)
=\frac{1}{\pi}\int_0^\infty e^{\lambda(\cos t-1)-\sigma^2t^2/2}\cos(\lambda\sin t+t(\mu-z))\,\mathrm dt.
$$
