Assume that $X$ and $Y$ follow the bivariate normal distribution with correlation coeffcient $\rho > 0$, zero means and scale parameters equal to one. I am looking for an elegant way to compute
$$E\left[ \Phi \left(X \right) \Phi \left(Y \right) \right]$$
I can see that the result is
$$E\left[ \Phi \left(X \right) \Phi \left(Y \right) \right] = \frac{1}{4} + \frac{\arcsin\left(\rho/2\right)} {\left(2\pi\right)}$$
but it's not obvious how one gets there. The CDFs are still uniformly distibuted RVs but the dependence complicates things. Had it not been for the dependence we would only have the $\frac{1}{4}$ factor on the RHS but now the shape of the distribution has to be taken into account.
I have experimented with the Law of Iterated Expectations but haven't gone very far. I don't think a trigonometric transformation would be helpful either. I would therefore appreciate it if someone could give some hints on how to approach this.
Thank you.