Evaluating $\int_0^\infty (t+a)^k e^{-t}\exp\left(-\frac{(t-\mu)^2}{2\sigma^2}\right)\,\mathrm dt$, $k\in\Bbb N_0$ 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)$.
 A: $$J_\ell=\int_0^\infty t^\ell e^{-t}\exp\left(-\frac{(t-\mu)^2}{2\sigma^2}\right)\,\,dt$$
Assuming $\Re(\mu )<\sigma ^2$, it seems that
$$J_\ell=2^{-\frac{\ell+1}{2}}\, e^{-\frac{\mu ^2}{2 \sigma ^2}} \,\sigma ^{\ell+1} \,\Gamma
   (\ell+1)\,\, U\left(\frac{\ell+1}{2},\frac{1}{2},\frac{\left(\mu -\sigma ^2\right)^2}{2
   \sigma ^2}\right)$$ where $U(a,b,c)$ is Tricomi's confluent hypergeometric function.
For integer values of $\ell$,  $U\left(\frac{\ell+1}{2},\frac{1}{2},t\right)$ is explicit and "quite" simple. We have
$$U\left(\frac{\ell+3}{2},\frac{1}{2},t\right)=\frac{2 \left((2 \ell+2 t+1) U\left(\frac{\ell+1}{2},\frac{1}{2},t\right)-2
   U\left(\frac{\ell-1}{2},\frac{1}{2},t\right)\right)}{(\ell+1) (\ell+2)}$$ (have a look at formula $13.3.7$ here)
A: A simpler integral is which is transformable into your is the following. The infinite radius of convergence makes the sum and integral switch possible: $$\int t^a e^{t^2+bt}dt=\int t^a\sum_{n=0}^\infty\frac{b^nt^n}{n!}\sum_{m=0}^\infty\frac{t^{2m}}{m!}dt= \sum_{n=0}^\infty\frac{b^n}{n!}\sum_{m=0}^\infty\frac1{m!}\int t^{a+2m+n}dt$$ to get: $$\int t^a e^{t^2+bt}dt =t^{a+1}\sum_{n=0}^\infty\sum_{m=0}^\infty\frac{b^nt^{2m+n}}{(2m+n+a+1)m!n!}$$
which works. Now use the lower incomplete gamma function $\gamma(a,z)$:
$$t^{a+1}\sum_{n=0}^\infty\sum_{m=0}^\infty\frac{b^nt^{2m+n}}{(2m+n+a+1)m!n!}=\left(-\frac1b\right)^{a+1}\sum_{n=0}^\infty\frac{\gamma(2n+a+1,-bt)}{b^{2n}n!}$$
which is simple enough, but no Wolfram Horn function nor Appell function match the series expansions.
A: 
Still working on it. Too big for a comment.

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