Gaussian Expectation of 2 Erfs 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?$$

Tanks for the hint. It was actually not that hard. I used 
$$\frac{u^2}{d}=x^2$$ leading to $$u = \sqrt{d} x$$ and $$du = \sqrt{d} dx.$$
Plugging-in yields
$$\int_{-\infty}^{\infty} exp{(-\frac{u^2}{d})}~ erf{(au)}~erf{(bu)}~du$$
$$=\int_{-\infty}^{\infty} exp{(-x^2)}~ erf{(a\sqrt{d} x)}~erf{(b\sqrt{d} x)}~\sqrt{d} dx$$
$$=\frac{2 \sqrt{d}}{\sqrt{\pi}} ~ atan\big(\frac{ab\sqrt{d}\sqrt{d}}{\sqrt{a^2d+b^2d+1}}\big)$$
$$=\frac{2 \sqrt{d}}{\sqrt{\pi}} ~ atan\big(\frac{ab\sqrt{d}}{\sqrt{a^2+b^2+\frac{1}{d}}}\big).$$
