I am interested in solving a definite double integral of the following form:
\begin{align} f(a,b) &= \int_0^\infty \exp\Big(\frac{-x^2}{2a}\Big)\int_{x}^{\infty} \exp\Big(\frac{-y^2}{2b}\Big) dydx\\ &= \int_0^\infty \exp\Big(\frac{-x^2}{2a}\Big)\text{erfc}\Big(\frac{x}{\sqrt{2b}}\Big)dx, \end{align}
for $a,b>0$, where erfc is the complementary error function. One potential way to go would be to use a power-series expansion (e.g., see answer by robjohn to this question or this paper), but I'm finding that a bit difficult to follow. I'm wondering if anybody has any ideas about ways to get an approximate answer. For now, I'm just trying to see if I can fit the numerical solution with a function of $a$ and $b$.
Consider $(X,Y)$ i.i.d. standard normal, then $$f(a,b)=2\pi\sqrt{ab}P(\sqrt bY>\sqrt aX>0)$$ The distribution of $(X,Y)$ is rotationally invariant hence, for every $c>0$, $$P(Y>cX>0)=\frac{\pi/2-\vartheta}{2\pi}$$ where $\vartheta$ in $(0,\pi/2)$ solves $$\tan\vartheta=c$$ Thus, $$f(a,b)=\sqrt{ab}(\pi/2-\arctan\sqrt{a/b})=\sqrt{ab}\arctan\sqrt{b/a}$$