The integral of a shifted error function with respect to the Gaussian measure admits a nice closed form expression:
$$ \int_{-\infty}^{\infty}\frac{1}{\sqrt{2\pi\sigma^{2}}}e^{-\frac{1}{2\sigma^{2}}(x-\mu)^{2}}{\rm erf}\left(ax+b\right) = {\rm erf}\left(\frac{a\mu+b}{\sqrt{1+2a^{2}\sigma^{2}}}\right) $$
I am after a similar expression for the squared error function: $$ \int_{-\infty}^{\infty}\frac{1}{\sqrt{2\pi\sigma^{2}}}e^{-\frac{1}{2\sigma^{2}}(x-\mu)^{2}}{\rm erf}\left(ax+b\right)^{2} =~ ??? $$ My first thought was integration by parts, which bring us to a Gaussian integral over $x{\rm erf}$, and for which another Math StackExchange has a rather complicated closed form expression in terms of Owen's T function.
This is a valid answer, but since I am not familiar with the properties of the Owen's T function, I was wondering if a simpler expression exists in this case, or at least for the simpler, particular case where $\mu=0$ and $\sigma=a=1$.
