I am struggling to approximate the following integral
$$\sqrt n\int_0^\infty \int_0^\infty (1 + n x^2)^{-1}(1 + y^2)^{-1} \Phi\left(\frac{a}{\sqrt{1 + b + x^2y^2}}\right) \text{d}x \text{d}y,$$ where $\Phi(u) = \frac{1}{2} \left\{1 + \text{erf}\left(\frac{u}{\sqrt2}\right)\right\}$ is the standard normal cumulative distribution function (and $\text{erf}$ is the error function). I also know that:
$n$ is a large natural number (so studying the integral in the asymptotic regime $n\rightarrow \infty$ can be sensible/useful);
$-4<a < -2$;
$b> 0$ (and it is close to $0$).
So far, I have been considering two directions to approximate this integral (but have been unsuccessful). I started by approximating the integral with respect to $x$:
using the series representation of the $\text{erf}$ function, but a) I feel that I would have to use many terms for the approximation to be accurate, and b) the integral is not much (?) simpler to compute when replacing $\Phi(\cdot)$ by the truncated series.
integrating by parts: an $\tan^{-1}$ term appears as well as the normal density function $\varphi(\cdot)$, and there I am stuck again...
Any idea to help me? In particular, I feel that using the fact that $n$ is large may help.
Thanks