Maybe this proof is too large, but you can deduce the result you want, sorry for that.

First of all:

- $\phi(x)= \frac{1}{\sqrt{2\pi}}e^{-\frac{x^2}{2}}$
- $\Phi(x)=\int_{-\infty}^{x}\phi(t)dt=\frac{1}{2}(1+erf(\frac{x}{\sqrt{2}}))$

To prepare the proof, start with this facts:

- $\int_{0}^{\infty}xe^{-x^2a}dx=\frac{1}{2a}$ (Proof by direct integral)
- $erf'(x)=\frac{2}{\sqrt{\pi}}e^{-x^2}$ (Proof by definition of $erf$ function)
- $$\int_0^{\infty}xe^{-\alpha^2x^2}erf(\beta x)dx=\frac{\beta}{2\alpha^2 \sqrt{\alpha^2+\beta^2}} $$
**Hint for 3.**
Integrating by parts, with
\begin{equation}u=erf(\beta x) \rightarrow du=\frac{2\beta}{\sqrt{\pi}}e^{-(x\beta)^2}dx\end{equation}
$$dv=xe^{-\alpha^2x^2} \rightarrow v=-\frac{e^{-x^2\alpha^2}}{2\alpha^2}$$
the result follows really easy.

Now, you have all the ingredients to calculate this monster:
$$T(b)=\int_{0}^{\infty}\Phi^2(bx)\phi(x)dx$$
Derivating $T(b)$ respect $b$, changing $\Phi(bx)$ by the definition in terms of $erf$ function, and using the fact 2, will lead you to:

$$T'(b)=...=\frac{1}{2\pi}\Big(\int_0^{\infty}{xe^{-x^2(b^2+1)/2}}dx + \int_0^{\infty}{x \cdot erf\big(\frac{bx}{\sqrt{2}}\big)e^{-x^2(b^2+1)/2}}dx \Big) =\frac{1}{2\pi}(A+B)$$

Compute A using the fact 1 to get:
$$A=\frac{1}{b^2+1}$$
and compute B using the fact 3 with $\alpha=\sqrt{\frac{b^2+1}{2}}$ and $\beta= \frac{b}{\sqrt{2}}$ to get:
$$B=\frac{b}{(b^2+1)\sqrt{2b^2+1}}.$$

It remains to integrate $T'(b)=\frac{A+B}{2\pi}$ respect $b$:
\begin{align}T(b)=&\frac{1}{2\pi}\Big(\int_{0}^{b}\frac{1}{t^2+1}dt+\int_{0}^{b}\frac{t}{(t^2+1)\sqrt{2t^2+1}}dt\Big)=\\
=&\frac{1}{2\pi}\Big(arctan(b)+arctan\big(\sqrt{2b^2+1}\big)\Big)+C. \hspace{0.5cm}\square
\end{align}

Mathematicaconfirms $T=7/24$ for $b=1$ but cannot solve the general case analytically. $\endgroup$ – David G. Stork Jul 18 '18 at 23:23