For a general case with $a,b$ we can derive the result of following integral, \begin{eqnarray} \int_{-\infty}^{\infty}\Phi(a+bx)^{2}\phi(x)dx\\&=&P\left(z_{1}\leq a+bx,z_{2}\leq a+bx\right)\\&=&P\left(z_{1}-bx\leq -a,z_{2}-bx\leq -a\right)\\&=& \mathcal{MVN}\left(x=\{-a,-a\},\mu=\{0,0\},\Sigma=\begin{bmatrix}b^{2}+1 & 1\\1& b^{2}+1 \end{bmatrix} \right) \end{eqnarray} In this case with $a=0, b=1$ we have, \begin{equation} \mathbb{E}(\Phi(X)^2)=\mathcal{MVN}\left(x=\{0,0\},\mu=\{0,0\},\Sigma=\begin{bmatrix}2 & 1\\1& 2 \end{bmatrix} \right) \end{equation}

After we take Jarle Tufto's result as given and having the fact that $\mathbb{V}(X)=1$, $\mathbb{V}(\Phi(X))=\mathbb{E}(\Phi(X)^2)-\mathbb{E}(\Phi(X))^2$ and $\mathbb{E}(\Phi(X))^2=\frac{1}{4}$, we then obtain final correlation formula \begin{eqnarray} \rho=\frac{\frac{1}{2\pi}}{\sqrt{1}\sqrt{\mathcal{MVN}\left(x=\{0,0\},\mu=\{0,0\},\Sigma=\begin{bmatrix}2 & 1\\1& 2 \end{bmatrix} \right)-\frac{1}{4}}} \end{eqnarray} A quick R-implementation shows,

sqrt(3/pi)
[1] 0.977205 

   (1/(2*sqrt(pi)))/sqrt(pmnorm(x = c(0,0),mean = rep(0.,2),varcov = matrix(c(2,1,1,2),ncol=2,byrow=T))-0.25)
[1] 0.977205

which coincides wolfies's result.

For a general case with $a,b$ we can derive the result of following integral, \begin{eqnarray} \int_{-\infty}^{\infty}\Phi(a+bx)^{2}\phi(x)dx\\&=&P\left(z_{1}\leq a+bx,z_{2}\leq a+bx\right)\\&=&P\left(z_{1}-bx\leq -a,z_{2}-bx\leq -a\right)\\&=& \mathcal{MVN}\left(x=\{-a,-a\},\mu=\{0,0\},\Sigma=\begin{bmatrix}b^{2}+1 & 1\\1& b^{2}+1 \end{bmatrix} \right) \end{eqnarray} In this case with $a=0, b=1$ we have, \begin{equation} \mathbb{E}(\Phi(X)^2)=\mathcal{MVN}\left(x=\{0,0\},\mu=\{0,0\},\Sigma=\begin{bmatrix}2 & 1\\1& 2 \end{bmatrix} \right) \end{equation}

After we take Jarle Tufto's result as given and having the fact that $\mathbb{V}(X)=1$, $\mathbb{V}(\Phi(X))=\mathbb{E}(\Phi(X)^2)-\mathbb{E}(\Phi(X))^2$ and $\mathbb{E}(\Phi(X))^2=\frac{1}{4}$, we then obtain final correlation formula \begin{eqnarray} \rho=\frac{\frac{1}{2\pi}}{\sqrt{1}\sqrt{\mathcal{MVN}\left(x=\{0,0\},\mu=\{0,0\},\Sigma=\begin{bmatrix}2 & 1\\1& 2 \end{bmatrix} \right)-\frac{1}{4}}} \end{eqnarray} A quick R-implementation shows,

    sqrt(3/pi)
    [1] 0.977205 

       (1/(2*sqrt(pi)))/sqrt(pmnorm(x = c(0,0), mean = rep(0.,2), 
            varcov = matrix(c(2,1,1,2), ncol=2, byrow=T))-0.25)
    [1] 0.977205

which coincides with wolfies's result.