Given that the joint density of $(X, Y)$ is $f(x, y) = \frac{1}{2\pi\sqrt{1 - \rho^2}}\exp\left(-\frac{x^2 + y^2 - 2\rho xy}{2(1 - \rho^2)}\right)$, the probability of interest is
\begin{align}
P(X > Y > 0) = \iint_{(x, y): x > y > 0}\frac{1}{2\pi\sqrt{1 - \rho^2}}\exp\left(-\frac{x^2 + y^2 - 2\rho xy}{2(1 - \rho^2)}\right)dxdy. \tag{1}\label{1}
\end{align}
A natural way of simplifying the integral $\eqref{1}$ is through the polar coordinates $x = \rho\cos\theta, y = r\sin\theta$, under which $\eqref{1}$ becomes
\begin{align}
P(X > Y > 0) = \frac{1}{2\pi\sqrt{1 - \rho^2}}\int_0^{\pi/4}\left[\int_0^\infty r\exp\left(-\frac{1 - \rho\sin(2\theta)}{2(1 - \rho^2)}r^2\right)dr\right]d\theta. \tag{2}\label{2}
\end{align}
It is easy to show that inner integral of $\eqref{2}$ evaluates to $\frac{1 - \rho^2}{1 - \rho\sin(2\theta)}$, whence it suffices to evaluate the integral
\begin{align}
I = \int_0^{\pi/4}\frac{1}{1 - \rho\sin(2\theta)}d\theta. \tag{3}\label{3}
\end{align}
One way to evaluate $\eqref{3}$ is considering the transformation $\theta = \arctan(u)$ (inspired by the tangent half-angle formula), which gives
\begin{align}
I = \int_0^1 \frac{1}{1 + u^2 - 2\rho u}du = \frac{1}{\sqrt{1 - \rho^2}}\arctan\left(\sqrt{\frac{1 + \rho}{1 - \rho}}\right).
\end{align}
Therefore,
\begin{align}
P(X > Y > 0) = \frac{1}{2\pi}\arctan\left(\sqrt{\frac{1 + \rho}{1 - \rho}}\right).
\end{align}
This agrees with the answer given in @Dilip Sarwate's comment (to verify, try matching derivatives of these two expressions).
An alternative approach makes use of the representation in your post. Assume $\rho \in (-1, 1)$ (the case $\rho = \pm 1$ is trivial), from the expression (I rewrote $X^i$ as $X_i$ to avoid confusion with exponents of $X$, and denote $\sqrt{\frac{1 - \rho}{1 + \rho}}$ by $\tau$ to save some typing)
$P(X > Y > 0) = P\left(X_2 > \tau X_1, X_1 > 0\right)$
that you have correctly derived and the condition that $X_1, X_2 \text{ i.i.d.} \sim N(0, 1)$, it follows that
\begin{align}
& P(X > Y > 0) = \int_0^\infty P\left(X_2 > \tau x\right)\phi(x)dx
= \int_0^\infty \left(1 - \Phi\left(\tau x\right)\right)\phi(x)dx \\
=& \frac{1}{2} - \int_0^\infty\int_{-\infty}^{\tau x}\phi(y)\phi(x) dydx \\
=& \frac{1}{2} - \frac{1}{2\pi}\int_0^\infty\int_{-\infty}^{\tau x}\exp\left(-\frac{1}{2}(x^2 + y^2)\right) dydx, \tag{4}\label{4}
\end{align}
where $\phi(x) = \frac{1}{\sqrt{2\pi}}e^{-x^2/2}$ and $\Phi(x) = \int_{-\infty}^x \phi(y)dy$.
It turns out that applying polar coordinates $x = r\cos\theta, y = r\sin\theta$ to $\eqref{4}$ would make the integration more straightforward (than to $\eqref{1}$). Indeed, by observing that the integration region $\{(x, y): x > 0, y < \tau x\}$ in the Cartesian coordinate system is the integration region $\{(r, \theta): r > 0, \cos\theta > 0, \sin\theta < \tau\cos\theta\} = \{(r, \theta): r > 0, \theta \in (0, \arctan\tau) \cup (\frac{3\pi}{2}, 2\pi)\}$ in the polar coordinate system, the integral in $\eqref{4}$ under polar coordinates becomes:
\begin{align}
\int_0^{\arctan\tau}\int_0^\infty re^{-\frac{1}{2}r^2}drd\theta + \int_{3\pi/2}^{2\pi}\int_0^\infty re^{-\frac{1}{2}r^2}drd\theta
= \arctan\tau + \frac{1}{2}\pi. \tag{5}\label{5}
\end{align}
Substituting $\eqref{5}$ back to $\eqref{4}$ yields
\begin{align}
P(X > Y > 0) = \frac{1}{4} - \frac{1}{2\pi}\arctan\tau
= \frac{1}{4} - \frac{1}{2\pi}\arctan\left(\sqrt{\frac{1 - \rho}{1 + \rho}}\right). \tag{6}\label{6}
\end{align}
Since $\arctan(x^{-1}) = \frac{\pi}{2} - \arctan(x)$ when $x > 0$, $\eqref{6}$ agrees with the answer found in the first approach.
