Vtyjkyuk

I am interested in solving a definite double integral of the following form:

beginalign
f(a,b) &= int_0^infty expBig(frac-x^22aBig)int_x^infty expBig(frac-y^22bBig) dydx\
&= int_0^infty expBig(frac-x^22aBig)texterfcBig(fracxsqrt2bBig)dx,
endalign

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$.

Assuming $a,b$ to be real and positive
$$I=int_0^infty e^-fracx^22 a, texterfcleft(fracxsqrt2b right),dx=sqrt2a int_0^infty e^-y^2 texterfcleft(sqrtfrac a b, yright),dy$$ Now, using series
$$texterfc(z)=1-frac2sqrtpisum_n=0^inftyfrac(-1)^n n!, (2n+1)z^2n+1$$ integrate termwise using the fact that
$$int_0^infty e^-z^2z^2n+1,dz=fracn!2$$ to get
$$I=sqrt2aleft(sqrtfrac pi 2-frac 1sqrt pisum_n=0^inftyfrac(-1)^n (2n+1)left(sqrtfrac a bright)^2n+1right)$$
$$I=sqrt2aleft(sqrtfrac pi 2-frac 1sqrt pitan ^-1left(sqrtfracabright)right)=sqrtfrac2api tan ^-1left(sqrtfracbaright)$$

Consider $(X,Y)$ i.i.d. standard normal, then $$f(a,b)=2pisqrtabP(sqrt bY>sqrt aX>0)$$ The distribution of $(X,Y)$ is rotationally invariant hence, for every $c>0$, $$P(Y>cX>0)=fracpi/2-vartheta2pi$$ where $vartheta$ in $(0,pi/2)$ solves $$tanvartheta=c$$ Thus, $$f(a,b)=sqrtab(pi/2-arctansqrta/b)=sqrtabarctansqrtb/a$$