Let $$ I_k=\int_0^\infty (t+a)^k e^{-t}\exp\left(-\frac{(t-\mu)^2}{2\sigma^2}\right)\,\mathrm dt, $$ with $k\in\Bbb N_0$ and $a>0$. Since $k$ is an integer we can expand the binomial to obtain $$ I_k=\sum_{\ell=0}^k\binom{k}{\ell}a^{k-\ell}\int_0^\infty t^\ell e^{-t}\exp\left(-\frac{(t-\mu)^2}{2\sigma^2}\right)\,\mathrm dt. $$ Expanding the quadratic in the Gaussian and combining all the exponential terms subsequently allows us to write a closed-form for $I_k$ that is a finite sum of parabolic cylinder function $D_\nu(z)$ with $$ D_\nu(z)=\frac{e^{-z^2/4}}{\Gamma(-\nu)}\int_0^\infty t^{-\nu-1} e^{-t^2/2-zt}\,\mathrm dt. $$
Can we write a closed-form for $I_k$ that does not involve a sum like this? Is there a special function, related to $D_\nu$, that admits an integral expression in the form of $I_k$? I would think that Meijer-G functions would be a potential candidate.
Edit:
I was asked for additional details/context. The origins of this problem are rooted in studying how photon noise passes through electro-optical image sensors. Without getting into too much detail, the model of the problem being studied leads to a random variable of the form $$ Y=\mathcal P(W)+R, $$ where $R\sim\mathcal N(0,\sigma_R^2)$ and $\mathcal P(W)$ is a compound Poisson random variable with random mean $W$, i.e. $\mathcal P(W)|W=w\sim\operatorname{Poisson}(w)$. In this problem, $W$ is truncated normal with lower bound $(a)$ and infinite upper bound. The density of $Y$ has the form $$ f_Y(y)=\sum_{k=0}^\infty \mathsf P(\mathcal P(W)=k)\phi(y-k,0,\sigma_R) $$ and the integral in question is needed to deduce $\mathsf P(\mathcal P(W)=k)$.