Computing the radioactive probability integral of a non-uniform system How do I integrate $e^{-y\alpha- x\beta- \gamma\sqrt{xy}}\,dy\,dx$, with $x$ and $y$ from $0$ to infinity, i.e., 
$$\int_0^\infty\int_0^\infty e^{-y\alpha- x\beta- \gamma\sqrt{xy}}\,dy\,dx\tag{1}$$ &
$$N\int_0^\infty\int_0^x e^{-y\alpha- x\beta- \gamma\sqrt{xy}}\,dy\,dx\tag{2}$$
where $N~=~(\text{equation } 1)^{-1}$

The change of variable $(2\beta x,2\alpha y)\to(x^2,y^2)$ shows that the integral to be computed is $$\iint_{x>0,y>0}e^{-\alpha y-\beta x-\gamma\sqrt{xy}}\,dxdy=\frac1{\alpha\beta}I\left(\frac\gamma{2\sqrt{\alpha\beta}}\right)$$ where, for every $|w|<1$, 
$$I(w)=\iint_{x>0,y>0}e^{-Q_w(x,y)/2}\,xydxdy$$ and $$Q_w(x,y)=x^2+y^2+2wxy$$
In particular, $I(w)=-J'(w)$, where $$J(w)=\iint_{x>0,y>0}e^{-Q_w(x,y)/2}\,dxdy$$
Since $|w|<1$, there exists $\vartheta$ in $(0,\pi)$ such that $$w=\cos\vartheta$$ The quadratic form $Q_w$ is positive, the inverse of the matrix of $Q_w$ being $$C=\begin{pmatrix}1&\cos\vartheta\\ \cos\vartheta&1\end{pmatrix}^{-1}=\frac1{\sin^2\vartheta}\begin{pmatrix}1&-\cos\vartheta\\-\cos\vartheta&1\end{pmatrix}$$ with $$|C|=\frac1{\sin^2\vartheta}$$
Thus, $$J(w))=\iint_{x>0,y>0}e^{-(x,y)^*C^{-1}(x,y)/2}dxdy=2\pi\sqrt{|C|}\,P(X>0,Y>0)$$ where $(X,Y)$ is centered normal with covariance $C$. Now, $$\sin\vartheta\, (X,Y)=(X_0,\sin\vartheta Y_0-\cos\vartheta X_0)$$ where $(X_0,Y_0)$ is standard normal hence $$[X>0,Y>0]=[X_0>0,Y_0>\cot\vartheta X_0]$$ The distribution of $(X_0,Y_0)$ is rotationally invariant hence the probability of this event is the angle of the sector $x>0$, $y>\cot\vartheta x$, divided by $2\pi$. This sector is limited by the angles $\frac\pi2-\vartheta$ and $\frac\pi2$ hence
$$P(X>0,Y>0)=\frac\vartheta{2\pi}$$ and $$J(w)=\frac\vartheta{\sin\vartheta}$$
Differentiating this with respect to $\vartheta$ and then $\vartheta$ with respect to $w$, one gets $$I(w)=\frac1{\sin^2\vartheta}\left(1-\vartheta\cot\vartheta\right)$$ Finally, for every $|w|<1$, $$I(w)=\frac1{1-w^2}\left(1-\frac{w}{\sqrt{1-w^2}}\arccos w\right)$$ In particular, $I(-1)$ is infinite, as was to be expected, and one can compute $$I(1)=\frac13$$
One could even deduce $I(w)$ for $w>1$ from this, by analytic continuation... but let us stop here.
