By a suitable transformation we reduce this to the case where $X$ and $Y$ have a bivariate normal distribution $\text{BN}(0,0,1,1,\rho)$ with joint density function $f(x,y;\rho)$.

We have

$$\tag{*}E(\Phi(X) \Phi(Y)) = \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} \left(\int_{-\infty}^x \phi(u)\, du\right) \left( \int_{-\infty}^y \phi(v)\, dv\right)f(x,y;\rho)\, dx \, dy,$$

where $\phi$ is the standard normal density function.

The key to solving this problem is to observe that (*) is equivalent to the joint probability that $U \leqslant X$ and $V \leqslant Y$ where $U$ and $V$ are standard normally distributed random variables that are uncorrelated with each other and where each is uncorrelated with $X$ and $Y$:

$$E(\Phi(X)\Phi(Y)) = P(U \leqslant X, V\leqslant Y) = P(X-U \geqslant 0, Y-V \geqslant 0).$$

Now $X-U$ and $Y-V$ are both normally distributed with mean $0$, standard deviation $\sqrt{2}$, and with correlation

$$\text{corr}(X-U,Y-V) = \frac{E((X-U)(Y-V))}{\sqrt{\text{var}(X-U)}\sqrt{\text{var}(Y-V)}} \\= \frac{E(XY)-E(XV) - E(YU) + E(UV)}{\sqrt{2}\sqrt{2}} \\ = \frac{\rho}{2}$$

It is shown here that if $X_1$ and $X_2$ have a standard bivariate normal distribution with correlation $\rho$, then

$$P(X_1 \geqslant 0, X_2 \geqslant 0) = \frac{1}{4} + \frac{\arcsin \rho}{2 \pi}.$$

Thus,

$$E(\Phi(X)\Phi(Y)) = \frac{1}{4} + \frac{1}{2\pi} \arcsin \frac{\rho}{2}.$$

Some further manipulation yields

$$\text{corr}(\Phi(X), \Phi(Y)) = \frac{6}{\pi} \arcsin \frac{\rho}{2}$$