I am considering the integral $\mathbb{E}[Z\Phi(aZ+b)^n]$, where $Z\sim\mathcal{N}(0,1)$, $\Phi$ is the cdf of a standard normal distribution, $n\in \mathbb{N}$, $a,b\in \mathbb{R}$. In the following, I let also $\phi$ denote the pdf of a standard normal distribution.
One approach, with inspiration from my previous question, reduced the integral to $$ an\int_{\mathbb{R}}\Phi(az+b)^{n-1}\phi(az+b) \phi(z)dz$$ If $a=1, b=0$, the two density functions could simplify, and after a simple scalation, I would be able to use the result from another of my recent questions.
If $n=1$, the result is on Wikipedia. However, here it is without proof, and it is also without proof in its reference Owen (1980), so I cannot get inspiration from there.
Another possible method would be to use a similar approach to the one in the previously mentioned question, but that only changed the $\Phi(aZ+b)^n$ to a multivariate cdf without seemingly helping much.
I hope some of you can help me out. Thank you very much!
You can show after some algebra that
$$\phi(z)\phi(az+b)=\frac{e^{-b^2/2(a^2+1)}}{\sqrt{2\pi}}\phi\left(\sqrt{a^2+1}~z+\frac{ab}{\sqrt{a^2+1}}\right)$$
Then after a change of variables the integral reduces to the expression
$$\mathbb{E}(Z\Phi(aZ+b)^n)=n\frac{a}{\sqrt{a^2+1}}\frac{\exp(-b^2/2(a^2+1))}{\sqrt{2\pi}}\mathbb{E}(\Phi(a'Z+b')^{n-1})$$
where $a':=\frac{a}{\sqrt{a^2+1}}~~,~~b':=\frac{b}{a^2+1}$
and this integral has been calculated in a linked question.