For the Standard Gaussian Expectation
$$\int_{-\infty}^{\infty} exp{(-u^2)}~ erf{(au)}~erf{(bu)}~du,$$
where $erf(x)$ is the cumulative densitiy function of a normal, I have found investigate the formulae
$$\frac{2}{\sqrt{\pi}} ~ atan\big(\frac{ab}{\sqrt{a^2+b^2+1}}\big)$$ (UPDATE: square-root was missing!)
described here. Can anybody confirm this result?
Anyway, I was wondering, how the formulae changes when instead a scaled Gaussian $exp{(-\frac{u^2}{d})}$ is used instead. That is, how can I compute
$$\int_{-\infty}^{\infty} exp{(-\frac{u^2}{d})}~ erf{(au)}~erf{(bu)}~du?$$