Integral of product of two error complementary functions (erfc) Could you please help me  to show that the integral
$$
\int_0^{\infty} \mathrm{erfc}(ax) \, \mathrm{erfc}(bx)\, \mathrm{d}x
$$
is equal to 
$$
\frac{1}{ab\sqrt{\pi}} (a+b-\sqrt{a^2+b^2}),
$$
where 
$$
\mathrm{erfc}(y)=\frac{2}{\sqrt{\pi}} \int_y^{\infty} \exp(-t^2)\, \mathrm{d} t.
$$
I have tried to expand the integral as
$$
\frac{4}{\pi }\int_0^{\infty} \int_{ax}^{\infty} \int_{bx}^{\infty} \exp(-t^2 -s^2) \, \mathrm{d}s \, \mathrm{d}t \, \mathrm{d} x
$$
but I could not come up with the right change of variables. Any ideas on how to proceed? 
Thank you in advance!
 A: 
Use the probabilistic way, which shows this is a problem of geometry in the plane...

That is, first note that, for every $x$, $$\mathrm{erfc}(x)=2P(X\geqslant \sqrt2 x)$$ where $X$ is standard normal hence, considering $(X,Y)$ i.i.d. standard normal and assuming that $a$ and $b$ are positive, one sees that the integral to be computed is $$I=4\int_0^\infty P(A_x)dx$$ where $$A_x:=[X\geqslant ax\sqrt2,Y\geqslant bx\sqrt2]$$ Now, $$(X,Y)=(R\cos\Theta,R\sin\Theta)$$ where $(R,\Theta)$ is independent, $\Theta$ is uniform on $(0,2\pi)$ and $R\geqslant0$ is such that $$P(R\geqslant r)=e^{-r^2/2}$$ for every $r\geqslant0$, hence $A_x\subset A_0=[\Theta\in(0,\pi/2)]$ and, on the event $A_0$, $$P(A_x\mid\Theta)=P(R\cos\Theta\geqslant ax\sqrt2,R\sin\Theta\geqslant bx\sqrt2\mid\Theta)=e^{-x^2u(\Theta)^2}$$ with $$u(\theta):=\max\left(\frac{a}{\cos\theta},\frac{b}{\sin\theta}\right)$$ Thus, still on $A_0$, $$\int_0^\infty P(A_x\mid\Theta)dx=\int_0^\infty e^{-x^2u^2(\Theta)}dx=\frac{\sqrt{\pi}}{2u(\Theta)}$$ This proves that $$\sqrt\pi I=4\sqrt\pi E\left(\int_0^\infty P(A_x\mid\Theta)dx\right)=2\pi E\left(\frac1{u(\Theta)}\mathbf 1_{A_0}\right)$$ that is, $$\sqrt\pi I=2\pi E\left(\min\left(\frac{\cos\Theta}a,\frac{\sin\Theta}b\right)\mathbf 1_{A_0}\right)$$ which is, using the distribution of $\Theta$, $$\sqrt\pi I=\int_0^{\pi/2}\min\left(\frac{\cos\theta}a,\frac{\sin\theta}b\right)d\theta=\frac1b\int_0^{\vartheta(a,b)}\sin\theta d\theta+\frac1a\int_{\vartheta(a,b)}^{\pi/2}\cos\theta d\theta$$ where the angle $\vartheta(a,b)$ in $(0,\pi/2)$ is uniquely defined by the condition that $$\tan\vartheta(a,b)=b/a$$ hence $$\sqrt\pi I=\frac{1-\cos\vartheta(a,b)}b+\frac{1-\sin\vartheta(a,b)}a$$ Finally, $$\cos\vartheta(a,b)=\frac{a}{\sqrt{a^2+b^2}}\qquad\sin\vartheta(a,b)=\frac{b}{\sqrt{a^2+b^2}}$$ hence $$I=\frac1{\sqrt{\pi}ab}\left(a-\frac{a^2}{\sqrt{a^2+b^2}}+b-\frac{b^2}{\sqrt{a^2+b^2}}\right)=\frac1{\sqrt{\pi}ab}\left(a+b-\sqrt{a^2+b^2}\right)$$
A: At first, we may remove the useless extra parameter, then study
$$ I(a) = \int_{0}^{+\infty}\text{erfc}(ax)\,\text{erfc}(x)\,dx $$
through differentiation under the integral sign. Assuming $a>0$ we have
$$ I'(a) = \frac{1}{\sqrt{\pi}}\int_{0}^{+\infty}-2x e^{-a^2 x^2}\text{erfc}(x)\,dx \stackrel{pp}{=}\frac{1}{a^2\sqrt{\pi}}-\frac{2}{\pi a^2}\int_{0}^{+\infty}e^{-(a^2+1)x^2}\,dx$$
from which:
$$ I'(a) = \frac{1}{\sqrt{\pi}}\left(\frac{1}{a^2}-\frac{1}{a^2\sqrt{1+a^2}}\right) $$
and:
$$ I(a) = \frac{1}{\sqrt{\pi}}\left(\frac{\sqrt{1+a^2}-1}{a}\right). $$
Through a substitution, if $a,b>0$ we get:
$$ \int_{0}^{+\infty}\text{erfc}(ax)\,\text{erfc}(bx)\,dx = \color{red}{\frac{a+b-\sqrt{a^2+b^2}}{ab\sqrt{\pi}}} $$
as wanted.
A: A multiple integral approach. Since
$$
\int_{0}^{\infty}x^2 e^{-\alpha x^2}\mathrm{d}x \overset{\alpha x^2 \to x}{=} \frac{1}{2\alpha^{\frac{3}{2}} }\int_{0}^{\infty} x^{\frac{3}{2}-1}e^{-x} \mathrm{d}x = \frac{1}{2\alpha^{\frac{3}{2}} } \Gamma\left(\frac{3}{2}\right)  = \frac{\sqrt{\pi}}{4\alpha^{\frac{3}{2}} } 
$$
We get
\begin{align}
\int_0^\infty \operatorname{erfc}(ax) \operatorname{erfc}(bx)\, \mathrm{d}x & =\int_0^\infty \left( \frac{2}{\sqrt{\pi}}\int_1^\infty ax e^{-a^2x^2z^2} \, \mathrm{d}z\right)\left( \frac{2}{\sqrt{\pi}}\int_1^\infty bx e^{-b^2x^2y^2}\, \mathrm{d}y\right) \mathrm{d}x\\
& =\frac{4ab}{\pi} \int_1^\infty \int_1^\infty \int_0^\infty x^2 e^{-(a^2z^2 +b^2y^2)x^2} \, \mathrm{d}x\, \mathrm{d}z\, \mathrm{d}y\\
& =\frac{ab}{\sqrt{\pi}} \int_1^\infty \int_1^\infty \frac{1}{\left(a^2z^2 +b^2y^2 \right)^{\frac32}} \, \mathrm{d}z\, \mathrm{d}y\\
& \overset{z\to \tan(\theta)by/a}{=}\frac{ab}{\sqrt{\pi}} \int_1^\infty \frac{1}{ab^2y^2} - \frac{1}{y^2b^2\sqrt{y^2b^2+a^2}} \, \mathrm{d}y\\
& =\frac{ab}{\sqrt{\pi}}\left(\frac{1}{ab^2} - \frac{1}{b^2}\int_1^\infty  \frac{1}{y^2\sqrt{b^2y^2+a^2}} \, \mathrm{d}y\right)\\
& \overset{y\to\tan(\varphi)a/b }{=}\frac{a}{b\sqrt{\pi}}\left(\frac{1}{a} - \frac{\sqrt{a^2+b^2} - b}{a^2}  \right)\\
& = \frac{a+b - \sqrt{a^2+b^2}}{ab\sqrt{\pi}}
\end{align}
