Can we evaluate $\int_0^\infty x^k\frac{ae^{bx}}{\left(1+ae^{bx}\right)^2}e^{-\frac{x^2}{2}}dx$ This sequence of integrals come from the following expectations $$E\left(X^kf(X)\right),\quad f(x;a,b)=\frac{ae^{bx}}{\left(1+e^{bx}\right)^2},\quad a>0,~b>0,~k=0,1,2,\cdots,$$
where $X\sim N(0,1)$. So we can express them as a sequence of integrals $$\int_{-\infty}^\infty x^k\frac{ae^{bx}}{\left(1+ae^{bx}\right)^2}e^{-\frac{x^2}{2}}dx,\quad k=0,1,2,\cdots.$$
Since $f(-x;a,b)=f(x;\frac{1}{a},b)$, which have some similarities, then we only consider $$\int_0^\infty x^k\frac{ae^{bx}}{\left(1+ae^{bx}\right)^2}e^{-\frac{x^2}{2}}dx.$$ I have found that $$\left(\frac{1}{b}\cdot\frac{ae^{bx}}{1+ae^{bx}}\right)'=\frac{ae^{bx}}{\left(1+ae^{bx}\right)^2},$$ so, maybe we can use integration by parts (here we take $k\geq 1$)
\begin{align}
\int_0^\infty x^k\frac{ae^{bx}}{\left(1+ae^{bx}\right)^2}e^{-\frac{x^2}{2}}dx&=\frac{1}{b}\int_0^\infty x^ke^{-\frac{x^2}{2}}d\frac{ae^{bx}}{1+ae^{bx}}\\
&=\frac{1}{b}x^ke^{-\frac{x^2}{2}}\frac{ae^{bx}}{1+ae^{bx}}\Big\vert_0^\infty-\frac{1}{b}\int_0^\infty \frac{ae^{bx}}{1+ae^{bx}}dx^ke^{-\frac{x^2}{2}}\\
&=\frac{1}{b}\int_0^\infty x^{k+1}\frac{ae^{bx}}{1+ae^{bx}}e^{-\frac{x^2}{2}}dx-\frac{k}{b}\int_0^\infty x^{k-1}\frac{ae^{bx}}{1+ae^{bx}}e^{-\frac{x^2}{2}}dx
\end{align}
Therefore, it suffices to calculate $$\int_0^\infty x^k\frac{ae^{bx}}{1+ae^{bx}}e^{-\frac{x^2}{2}}dx,\quad k=0,1,2\cdots,$$
and I am stuck on this.
 A: The farthest I could go
$$x^k\frac{ae^{bx}}{1+ae^{bx}}e^{-\frac{x^2}{2}}=\sum_{n=0}^\infty (-1)^n\, x^k\,e^{-\frac{x^2}{2}}\left(\frac{e^{-b x}}{a}\right)^n$$
Defining
$$J_{n,k}=\int_0^\infty x^k\,e^{-\frac{x^2}{2}}e^{-nb x}\,dx$$ Using Kummer confluent hypergeometric functions, they write (with $t=bn$)
$$J_{n,k}=2^{\frac{k-1}{2}} \Gamma
   \left(\frac{k+1}{2}\right) \,
   _1F_1\left(\frac{k+1}{2};\frac{1}{2};\frac{t^2
   }{2}\right)- 2^{\frac k2} t\, \Gamma
   \left(\frac{k+2}{2}\right) \,
   _1F_1\left(\frac{k+2}{2};\frac{3}{2};\frac{t^2}
  {2}\right)$$
They write
$$J_{n,k}=(-1)^k\,\sqrt{\frac{\pi }{2}} e^{\frac{t^2 }{2}}
   \text{erfc}\left(\frac{t}{\sqrt{2}}\right)\,P_k(t)+(-1)^{k+1}\,Q_{k}(t)$$ and the first polynomials are respectively
$$\left(
\begin{array}{cc}
 k & P_k(t) \\
 0 & 1 \\
 1 & t \\
 2 & t^2+1 \\
 3 & t^3+3 t \\
 4 & t^4+6 t^2+3 \\
 5 & t^5+10 t^3+15 t \\
 6 & t^6+15 t^4+45 t^2+15 \\
 7 & t^7+21 t^5+105 t^3+105 t \\
 8 & t^8+28 t^6+210 t^4+420 t^2+105 \\
\end{array}
\right)$$
and
$$\left(
\begin{array}{cc}
 k & Q_k(t) \\
 0 & 0 \\
 1 & 1 \\
 2 & t \\
 3 & t^2+2 \\
 4 & t^3+5 t \\
 5 & t^4+9 t^2+8 \\
 6 & t^5+14 t^3+33 t \\
 7 & t^6+20 t^4+87 t^2+48 \\
 8 & t^7+27 t^5+185 t^3+279 t \\
\end{array}
\right)$$
These show interesting patterns which would be intresting to explore further.
The problem remains with the infinite summations.
