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)$.
Well, we are trying to solve the following integral:
$$\mathcal{I}_\text{k}\left(\alpha,\mu,\sigma\right):=\int\limits_{0}^{\infty}\left(x+\alpha\right)^\text{k}\exp\left(-x\right)\exp\left(-\frac{1}{2}\cdot\left(\frac{x-\mu}{\sigma}\right)^2\right)\space\text{d}x\tag1$$
Using the binomial expansion we can write:
$$\mathcal{I}_\text{k}\left(\alpha,\mu,\sigma\right)=\sum_{\text{n}\space=\space0}^\text{k}\binom{\text{k}}{\text{n}}\alpha^{\text{k}-\text{n}}\int\limits_{0}^{\infty}x^\text{n}\exp\left(-x\right)\exp\left(-\frac{1}{2}\cdot\left(\frac{x-\mu}{\sigma}\right)^2\right)\space\text{d}x\tag2$$
Now, we can see that using Laplace transforms:
$$\mathcal{I}_\text{k}\left(\alpha,\mu,\sigma\right)=\lim_{\text{s}\space\to\space1}\sum_{\text{n}\space=\space0}^\text{k}\binom{\text{k}}{\text{n}}\alpha^{\text{k}-\text{n}}\int\limits_{0}^{\infty}x^\text{n}\exp\left(-\text{s}x\right)\exp\left(-\frac{1}{2}\cdot\left(\frac{x-\mu}{\sigma}\right)^2\right)\space\text{d}x\tag3$$
Using the properties of the Laplace transform, we can write:
$$\mathcal{I}_\text{k}\left(\alpha,\mu,\sigma\right)=\lim_{\text{s}\space\to\space1}\sum_{\text{n}\space=\space0}^\text{k}\binom{\text{k}}{\text{n}}\left(-1\right)^\text{n}\alpha^{\text{k}-\text{n}}\cdot\frac{\partial^\text{n}}{\partial\text{s}^\text{n}}\left(\mathscr{L}_x\left[\exp\left(-\frac{1}{2}\cdot\left(\frac{x-\mu}{\sigma}\right)^2\right)\right]_{\left(\text{s}\right)}\right)\tag4$$
Now, we can see that:
$$\exp\left(-\frac{1}{2}\cdot\left(\frac{x-\mu}{\sigma}\right)^2\right)=\exp\left(-\frac{1}{2}\cdot\left(\frac{x}{\sigma}\right)^2\right)\exp\left(\frac{x\mu}{\sigma^2}\right)\exp\left(-\frac{1}{2}\cdot\left(\frac{\mu}{\sigma}\right)^2\right)\tag5$$
So, we get:
$$\displaystyle\mathcal{I}_\text{k}\left(\alpha,\mu,\sigma\right)=\exp\left(-\frac{1}{2}\cdot\left(\frac{\mu}{\sigma}\right)^2\right)\lim_{\text{s}\space\to\space1}\sum_{\text{n}\space=\space0}^\text{k}\binom{\text{k}}{\text{n}}\left(-1\right)^\text{n}\alpha^{\text{k}-\text{n}}\cdot\frac{\partial^\text{n}}{\partial\text{s}^\text{n}}\left(\mathscr{L}_x\left[\exp\left(-\frac{1}{2}\cdot\left(\frac{x}{\sigma}\right)^2\right)\exp\left(\frac{x\mu}{\sigma^2}\right)\right]_{\left(\text{s}\right)}\right)\tag6$$