Double integral involving the normal CDF I need to compute (or best approximate?) the following integral
$$\int_0^\infty \int_0^\infty (1 + \alpha u)^{-1}(1 + v)^{-1} \Phi\left(\frac{\beta}{\sqrt{\gamma + uv}}\right) \text{d}u \text{d}v,\qquad \gamma > 0, \quad \beta \in \mathbb{R}, \quad \alpha \in \mathbb{N}^*,$$
where $\Phi(\cdot)$ is the standard normal cumulative distribution function and $\alpha$ is typically a large natural number. What strategy would you advise for this?
I don't know if this is helpful but I can also reformulate my problem (up to multiplicative constants) as computing the following expectation:
$$\text{E}\left\{ \Phi\left(\frac{\beta}{\sqrt{\gamma + X^2Y^2}}\right)\right\},\qquad \gamma > 0, \quad \beta \in \mathbb{R},$$
where 
$$X \sim \text{Half-Cauchy}(0,1), \qquad Y \sim \text{Half-Cauchy}(0,\alpha^{-1/2}), \quad \alpha \in \mathbb{N}^* \, (\text{large}).$$
EDIT: I made a mistake in the first integral, it should be:
$$\int_0^\infty \int_0^\infty (1 + \alpha x^2)^{-1}(1 + y^2)^{-1} \Phi\left(\frac{\beta}{\sqrt{\gamma + x^2y^2}}\right) \text{d}x \text{d}y,\quad \gamma > 0, \beta \in \mathbb{R}, \alpha \in \mathbb{N}^*\, \text{(large)},$$
which now corresponds to the above expectation.
EDIT 2: Another way to solve my problem would be to characterise the difference (hopefully small?)
$$\text{E}\left\{ \Phi\left(\frac{\beta}{\sqrt{\gamma + X^2Y^2}}\right)\right\} - \Phi\left(\frac{\beta}{\sqrt{\gamma}}\right),\qquad \gamma > 0, \quad \beta \in \mathbb{R},$$
would it be sensible to approximate the difference
$$\Phi\left(\frac{\beta}{\sqrt{\gamma + X^2Y^2}}\right) - \Phi\left(\frac{\beta}{\sqrt{\gamma}}\right),$$
with a Taylor series in zero or is this unreasonable given that the half-Cauchy distributions are heavy-tailed? 
 A: Under the conditions that $\alpha, \gamma, u,$ and $v$ are all positive,
$$\beta/\sqrt{\gamma + uv} \ge \min(\beta/\sqrt{\gamma}, 0) = \delta \gt -\infty.$$
Therefore, because $\Phi$ is a CDF for a distribution supported on $(-\infty,\infty),$ $$\Phi\left(\frac{\beta}{\sqrt{\gamma + uv}}\right)\ge \Phi(\delta) = \epsilon \gt 0.$$
Consequently
$$\eqalign{
&\int_0^\infty \int_0^\infty (1+\alpha u)^{-1}(1+v)^{-1} \Phi\left(\frac{\beta}{\sqrt{\gamma + uv}}\right) \mathrm{d}u\mathrm{d}v \\
&\ge 
\epsilon\int_0^\infty \int_0^\infty (1+\alpha u)^{-1}(1+v)^{-1}  \mathrm{d}u\mathrm{d}v \\
&= \lim_{M\to\infty}\lim_{N\to\infty}\epsilon\int_0^M (1+\alpha u)^{-1}\mathrm{d}u\int_0^N (1+v)^{-1} \mathrm{d}v \\
&=\frac{\epsilon}{\alpha} \lim_{M\to\infty}\lim_{N\to\infty} \log(1 + M\alpha)\log(1 + N),
}$$ 
which diverges to $+\infty.$
A: If approximation is OK, you could simulate the expected value through an average. A law of large numbers implies that the average converges to the expected value - provided the latter exists, which is however not obvious to me. I did simulate a few runs, though, and got very similar results each time, which is not indicative of an issue with heavy tails.
library(LaplacesDemon)

draws <- 1e6
alpha <- 3
gamma <- 1
beta <- 1
x <- rhalfcauchy(draws, scale=1)
y <- rhalfcauchy(draws, scale=1/sqrt(alpha))

mean(pnorm(beta/sqrt(gamma+x^2*y^2)))

