Closed form for $\int_{0}^{\infty }\!{\rm erf} \left(cx\right) \left( {\rm erf} \left(x \right) \right) ^{2}{{\rm e}^{-{x}^{2}}}\,{\rm d}x$ I encountered this integral in my calculations:
$$\int_{0}^{\infty }\!{\rm erf} \left(cx\right) \left( {\rm erf} \left(x
\right) \right) ^{2}{{\rm e}^{-{x}^{2}}}\,{\rm d}x$$
where $c>0$ and $c\in \mathbb{R}$
but could not find a closed-form representation for it.
I also tried to find possible closed forms using Inverse Symbolic Calculator and WolframAlpha but they did not find anything.
I was looking in the book "Integrals ans Series Volume 2-Prudnikov-Brychkov-Marychev" but did not find a similar formula. 
I am not sure it exists, but if it exists I want to know it. Closed forms are easier to manipulate, sometimes closed forms of different integrals or sums contain terms that cancel each other etc
Could you please help me to find a closed form (even using non-elementary special functions), if it exists?
 A: I think it is more practical to consider the parametric integral
$$ I(a,b,c) = \int_{0}^{+\infty}\!\!\!\text{erf}(ax)\,\text{erf}(bx)\,\text{erf}(cx)\,e^{-x^2}dx \tag{1}$$
and notice that:
$$ \frac{\partial^3 I}{\partial a\,\partial b\,\partial c}=\frac{8}{\pi^{3/2}}\int_{0}^{+\infty}x^3 e^{-(1+a^2+b^2+c^2)x^2}\,dx=\frac{4}{\pi^{3/2}(1+a^2+b^2+c^2)^2}\tag{2}$$
Integrating over $(b,c)\in(0,1)^2$ then with respect to $a$ is now tedious but doable with the support of a CAS. 
Obviously, we would be happier to integrate over $0\leq b^2+c^2\leq 2$ or $0\leq b^2+c^2\leq 1$ for first: this change of integration domain leads to decent bounds for $I(a,b,c)$. Anyway, the problem is equivalent to integrating a Cauchy-like distribution over $R=[0,1]\times[0,1]\times[0,a]$: a probabilistic approach through characteristic functions might be interesting, too.
A: Partial work. Consider $$I\left(c\right)=\int_{0}^{\infty}\textrm{erf}\left(cx\right)\textrm{erf}^{2}\left(x\right)e^{-x^{2}}dx.
 $$ Then we have, integrating by parts, that $$I\left(c\right)=\frac{\sqrt{\pi}}{6}-\frac{c}{3}\int_{0}^{\infty}\textrm{erf}^{3}\left(x\right)e^{-c^{2}x^{2}}dx=\frac{\sqrt{\pi}}{6}-\frac{c}{3}\int_{0}^{\infty}\left(1-\textrm{erfc}\left(x\right)\right)^{3}e^{-c^{2}x^{2}}dx
 $$ $$=\frac{\sqrt{\pi}}{6}-\frac{c}{3}\int_{0}^{\infty}e^{-c^{2}x^{2}}dx+\frac{c}{3}\int_{0}^{\infty}\textrm{erfc}^{3}\left(x\right)e^{-c^{2}x^{2}}dx
 $$ $$-c\int_{0}^{\infty}\textrm{erfc}^{2}\left(x\right)e^{-c^{2}x^{2}}dx+c\int_{0}^{\infty}\textrm{erfc}\left(x\right)e^{-c^{2}x^{2}}dx
 $$ $$=\frac{\sqrt{\pi}}{6}-\frac{c}{3}H_{0,0}\left(c^{2}\right)+\frac{c}{3}H_{0,3}\left(c^{2}\right)-cH_{0,2}\left(c^{2}\right)+cH_{0,1}\left(c^{2}\right),
 $$ say. I used this particular notation because now we may use these results at page $13$ and get $$H_{0,0}\left(c^{2}\right)=\frac{\sqrt{\pi}}{2c}
 $$ $$H_{0,1}\left(c^{2}\right)=\frac{\arctan\left(c\right)}{c\sqrt{\pi}}
 $$ $$H_{0,2}\left(c^{2}\right)=\frac{1}{c\sqrt{\pi}}\left(2\arctan\left(c\right)-\arccos\left(\frac{1}{1+c^{2}}\right)\right)
 $$ and, using the corollary $10.7$, $$H_{0,3}\left(c^{2}\right)=\frac{\sqrt{\pi}}{2c}-\frac{\arctan\left(1/c\right)}{c\sqrt{\pi}}+\frac{6}{c\pi}\int_{c^{2}}^{\infty}\frac{\arctan\left(\sqrt{t+2}\right)}{\sqrt{t+2}\sqrt{t}\left(t+1\right)}dt.
 $$ At this moment I don't know if it is possible to evaluate the last integral, I will work on it later. 
However it is quite simple to prove that $$\int_{c^{2}}^{\infty}\frac{1}{\sqrt{t+2}\sqrt{t}\left(t+1\right)}dx=\frac{\pi}{2}-2\arctan\left(\frac{c}{\sqrt{c^{2}+2}}\right)
 $$ so obviously $$\int_{c^{2}}^{\infty}\frac{\arctan\left(\sqrt{t+2}\right)}{\sqrt{t+2}\sqrt{t}\left(t+1\right)}dt\leq\frac{\pi}{4}-\pi\arctan\left(\frac{c}{\sqrt{c^{2}+2}}\right)
 $$ and $$\int_{c^{2}}^{\infty}\frac{\arctan\left(\sqrt{t+2}\right)}{\sqrt{t+2}\sqrt{t}\left(t+1\right)}dt\geq\left(\frac{\pi}{2}-2\arctan\left(\frac{c}{\sqrt{c^{2}+2}}\right)\right)\arctan\left(\sqrt{c^{2}+2}\right)
 $$ so it is not a closed form but it could help.
