We can first evaluate the probability $P[U > a]$ by noting $X^2 + Y^2 \overset{d}{=} \xi^2 + \eta^2 + 2\rho\xi\eta$, where $\xi, \eta$ i.i.d. $\sim N(0, 1)$. Therefore, for $a > 0$:
\begin{align}
& P[X^2 + Y^2 > a] = P[\xi^2 + \eta^2 + 2\rho\xi\eta > a] \\
=& \iint\limits_{[(x, y): x^2 + y^2 + 2\rho xy > a]}\frac{1}{2\pi}\exp\left(-\frac{1}{2}(x^2 + y^2)\right)dxdy. \tag{1}
\end{align}
To evaluate this double integral, apply the polar coordinates transformation $x = r\cos\theta, y = r\sin\theta$ with $r > 0, \theta \in [0, 2\pi)$. The integral $(1)$ then becomes
\begin{align}
\iint\limits_{[(r, \theta): (1 + \rho\sin(2\theta))r^2 > a]}\frac{1}{2\pi}r\exp\left(-\frac{1}{2}r^2\right)drd\theta. \tag{2}
\end{align}
Note the region $[(r, \theta): (1 + \rho\sin(2\theta))r^2 > a]$ is contained in the region $[(r, \theta): (1 + |\rho|)r^2 > a]$ and the integrand is positive, hence the integral $(2)$ is bounded above by
\begin{align}
\iint\limits_{[(r, \theta): (1 + |\rho|)r^2 > a]}\frac{1}{2\pi}r\exp\left(-\frac{1}{2}r^2\right)drd\theta
= \int_{\sqrt{\frac{a}{1 + |\rho|}}}^\infty e^{-r^2/2}rdr =
\exp\left(-\frac{a}{2(1 + |\rho|)}\right).
\end{align}
This completes the proof.
To see why $X^2 + Y^2 \overset{d}{=} \xi^2 + \eta^2 + 2\rho\xi\eta$, note that if $\begin{bmatrix} X \\ Y \end{bmatrix} \sim N_2(0, \Sigma)$, then $\begin{bmatrix} X \\ Y \end{bmatrix} \overset{d}{=} C\begin{bmatrix} \xi \\ \eta \end{bmatrix}$, where $\begin{bmatrix} \xi \\ \eta \end{bmatrix} \sim N_2(0, I_{(2)})$ and $C = \Sigma^{1/2}$ is the square root matrix of $\Sigma$. It then follows that
\begin{align}
X^2 + Y^2 &= \begin{bmatrix} X & Y \end{bmatrix}\begin{bmatrix} X \\ Y \end{bmatrix} \\
& \overset{d}{=} \begin{bmatrix} \xi & \eta \end{bmatrix}C'C\begin{bmatrix} \xi \\ \eta \end{bmatrix} = \begin{bmatrix} \xi & \eta \end{bmatrix}\Sigma\begin{bmatrix} \xi \\ \eta \end{bmatrix} \\
&= \xi^2 + \eta^2 + 2\rho\xi\eta.
\end{align}
