Analytic solution for integral of Gaussian and error functions

I am trying to find an analytic solution to the following definite integral:

$$\int_{-\infty}^\infty \mathrm{d}x~x^2 \Phi(a+bx)^2 \phi(x)$$

where $\Phi(a+bx) = \frac{1}{2}\left[1 + \mathrm{erf}\left(\frac{a+bx}{\sqrt{2}}\right)\right]$ and $\phi(x) = \frac{\exp\left[-\frac{x^2}{2}\right]}{\sqrt{2\pi}}$. I have attempted integrating by parts using some of the integrals on this table with no luck (though there may be some combinations I missed). Does anyone know if an analytic solution exists for this integral and what the solution might be?

I managed to come up with a solution using a combination of integration by parts and variable substitution. Here is the solution: $$J = \int_{-\infty}^\infty \mathrm{d}x~x^2 \Phi(a+bx)^2 \phi(x) = J_1+J_2$$ where $$J_1 = \frac{2b}{1+b^2}\phi\left(\frac{a}{\sqrt{1+b^2}}\right)\left(K_1 - K_2\right),$$ $$K_1 = \frac{b}{\sqrt{1+2b^2}} \phi\left(\frac{a}{\sqrt{(1+b^2)(1+2b^2)}}\right),$$ $$K_2 = \frac{ab}{\sqrt{1+b^2}} \Phi\left(\frac{a}{\sqrt{(1+b^2)(1+2b^2)}}\right),$$ $$J_2 = \Phi\left(\frac{a}{\sqrt{1+b^2}}\right) - 2\mathcal{T}\left(\frac{a}{\sqrt{1+b^2}}, \frac{1}{\sqrt{1+2b^2}}\right),$$ and $\mathcal{T}$ is Owen's T function.
Here is how I solved it. The following quantities will appear in some of the integrals: $$a' = \frac{ab}{\sqrt{1+b^2}}, \ b' = \sqrt{1+b^2}, \ a'' = \frac{ab'-a'b}{b'}, \text{ and } b'' = \frac{b}{b'}.$$ The goal at each step is to reduce each integral to one of the integrals with known solution from the table here.
(1) Integrate $J$ by parts with $u = \Phi(a+bx)^2 \Rightarrow \mathrm{d}u = 2b\phi(a+bx)\Phi(a+bx)~\mathrm{d}x$ and $\mathrm{d}v = x^2 \phi(x)~\mathrm{d}x \Rightarrow v = \Phi(x)-x\phi(x)$: $$J = 1 - 2b\int_{-\infty}^\infty \mathrm{d}x~\left[\Phi(x) - x\phi(x)\right]\phi(a+bx)\Phi(a+bx) = J_1+J_2$$ where $$J_1 = 2b\int_{-\infty}^\infty \mathrm{d}x~x\phi(x)\phi(a+bx)\Phi(a+bx)$$ and $$J_2 = 1-2b\int_{-\infty}^\infty \mathrm{d}x~\phi(a+bx)\Phi(a+bx)\Phi(x).$$ (2) Integrate $J_2$ by parts with $u = \Phi(x) \Rightarrow \mathrm{d}u = \phi(x)~\mathrm{d}x$ and $\mathrm{d}v = \phi(a+bx) \Phi(a+bx)~\mathrm{d}x \Rightarrow v = \frac{1}{2b}\Phi(a+bx)^2$: $$J_2 = \int_{-\infty}^\infty \mathrm{d}x~\Phi(a+bx)^2 \phi(x) = \Phi\left(\frac{a}{\sqrt{1+b^2}}\right) - 2\mathcal{T}\left(\frac{a}{\sqrt{1+b^2}}, \frac{1}{\sqrt{1+2b^2}}\right).$$
(3) With respect to $J_1$, the product of the Gaussians is rewritten as $\phi(x)\phi(a+bx) = \phi(a'+b'x)\phi\left(\frac{a'}{b}\right)$ and then the variable substitution $y = a' + b'x$ is employed: $$J_2 = \frac{2b}{b^{'2}} \phi\left(\frac{a'}{b}\right) \int_{-\infty}^\infty \mathrm{d}y~(y-a')\phi(y)\Phi(a''+b''y) = \frac{2b}{b^{'2}} \phi\left(\frac{a'}{b}\right) \left(K_1 - K_2\right)$$ where $$K_1 = \int_{-\infty}^\infty \mathrm{d}y~y\phi(y)\Phi(a''+b''y) = \frac{b''}{\sqrt{1+b^{''2}}} \phi\left(\frac{a''}{\sqrt{1+b^{''2}}}\right)$$ and $$K_2 = a' \int_{-\infty}^{\infty} \mathrm{d}y~\phi(y)\Phi(a''+b''y) = a' \Phi\left(\frac{a''}{\sqrt{1+b^{''2}}}\right)$$
(4) Now all terms have been integrated and all that is left is to substitute $J_1$, $J_2$, $K_1$, $K_2$, $a'$, $b'$, $a''$, and $b''$ back into the equation to recover $J$ in terms of $a$ and $b$ only.