Let us write $$\mathrm{erf}(x)=\frac{2}{\sqrt {\pi}}\int_0^x e^{-t^2}dt $$ for the usual Gauss error function. Given natural numbers $m,n,k$ I am interested in computing the integral $$\int_{-\infty}^{\infty} x^k\left(1+\mathrm{erf} (x)\right)^m\left(1-\mathrm{erf} (x)\right)^n e^{-x^2} dx. $$ The case $k=0$ is fairly straightforward and yields a result involving the Euler Beta function. The cases $k=1,2,... $ however seem to be a bit more difficult. I am particularly interested in the case $k=2$.
A related integral that I am equally interested in is $$\int_{-\infty}^{\infty}\left(1+\mathrm{erf} (x)\right)^m\left(1-\mathrm{erf} (x)\right)^n e^{-\alpha x^2} dx $$ with $\alpha>0$. It would be sufficient to obtain a closed form of the integral $$\int_{0}^{\infty}\left[\mathrm{erf} (x)\right]^m e^{-\alpha x^2} dx.$$ Here the case $\alpha=1$ admits a rather straightforward solution, while in general things seem to be a bit more involved.