Suppose, $X_1$ and $X_2$ are independent inverse Gaussian random variables with parameters $(\mu, \lambda_1)$ and $(\mu, \lambda_2)$, respectively. (See Wikipedia for the notation).
I am interested in the following probability (and whether or not it has a closed form solution):
$$ P(X_1 < X_2) = \int_0^\infty \int_0^y \exp \left(- \frac{\lambda_1(x-\mu)^2}{2\mu^2x} - \frac{\lambda_2(y-\mu)^2}{2\mu^2y} \right) \sqrt{\frac{\lambda_1 \lambda_2}{(2\pi)^2 x^3 y^3}} dx dy$$
$$= \int_0^\infty e^{-\frac{\lambda_2(y-\mu)^2}{2\mu^2y}}\left(\Phi \left(\sqrt{\frac{\lambda_1}{y}}\left(\frac{y}{\mu} -1\right) \right) + e^{2\lambda_1/\mu}\Phi \left(-\sqrt{\frac{\lambda_1}{y}}\left(\frac{y}{\mu} +1\right) \right) \right) \sqrt{\frac{\lambda_2}{2\pi y^3}}dy$$
where $\Phi$ is the cdf (Cumulative distribution function), of a standard normal random variable. I have no idea how to deal with this integral. The reason I am interested in its closed form is due to the fact that for the case $\mu = 0$, there is a very nice closed form solution relating to the hitting times of two independent Brownian motions (see here).
I have tried consulting this table for hints of how to deal with the integral of the erf function, against some other complicated term, to no avail. Any help would be massively appreciated!