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.
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.