How to evaluate $\int_0^\infty\operatorname{erfc}^n x\ \mathrm dx$? Let $\operatorname{erfc}x$ be the complementary error function.
I successfully evaluated these integrals:
$$\int_0^\infty\operatorname{erfc}x\ \mathrm dx=\frac1{\sqrt\pi}\tag1$$
$$\int_0^\infty\operatorname{erfc}^2x\ \mathrm dx=\frac{2-\sqrt2}{\sqrt\pi}\tag2$$
(Both $\operatorname{erfc}x$ and $\operatorname{erfc}^{2}x$ have primitive functions in terms of the error function.)
But I have problems with
$$\int_0^\infty\operatorname{erfc}^3x\ \mathrm dx\tag3$$
and a general case
$$\int_0^\infty\operatorname{erfc}^n x\ \mathrm dx.\tag4$$
Could you suggest an approach to evaluate them as well?
 A: Integrating by parts twice, we get
$$ \begin{align} \int_{0}^{\infty} \text{erfc}^{3}(x) \, dx &= x \,  \text{erfc}^{3}(x) \Bigg|^{\infty}_{0} + \frac{6} {\sqrt{\pi}} \int_{0}^{\infty}x \, \text{erfc}^{2}(x) e^{-x^{2}} \, dx  \\ &= \frac{6} {\sqrt{\pi}} \int_{0}^{\infty}x \, \text{erfc}^{2}(x) e^{-x^{2}} \, dx \\ &= - \frac{3 \, \text{erfc}(x) e^{-x^{2}}}{\sqrt{\pi}} \Bigg|_{0}^{\infty} - \frac{12}{\pi} \int_{0}^{\infty}\text{erfc}(x) e^{-2x^{2}} \, dx \\ &= \frac{3}{\sqrt{\pi}} - \frac{12}{\pi} \int_{0}^{\infty} \text{erfc}(x) e^{-2x^{2}} \, dx. \end{align}$$
And using the integral definition of the complementary error function, we get
$$ \begin{align} \int_{0}^{\infty} \text{erfc}(x) e^{-2x^{2}} \, dx &= \frac{2}{\sqrt{\pi}} \int_{0}^{\infty} \int_{1}^{\infty} x e^{-x^{2}t^{2}} e^{-2x^{2}} \, dt \, dx \\ &= \frac{2}{\sqrt{\pi}} \int_{1}^{\infty} \int_{0}^{\infty} x e^{-(2+t^{2})x^{2}} \, dx \, dt \\ &= \frac{1}{\sqrt{\pi}} \int_{1}^{\infty} \frac{1}{2+t^{2}} \, dt \\ &= \frac{\sqrt{2}}{2 \sqrt{\pi}} \int_{1 / \sqrt{2}}^{\infty} \frac{1}{1+u^{2}} \, du \\ &= \frac{\sqrt{2}}{2 \sqrt{\pi}} \left[ \frac{\pi}{2} - \arctan \left(\frac{1}{\sqrt{2}} \right) \right] \\ &= \frac{\sqrt{2}}{2 \sqrt{\pi}} \,  \arctan (\sqrt{2}). \end{align}$$
Therefore,
$$ \begin{align} \int_{0}^{\infty} \text{erfc}^{3}(x) \ dx &= \frac{3}{\sqrt{\pi}} - \frac{12}{\pi} \left( \frac{\sqrt{2}}{2 \sqrt{\pi}} \arctan \sqrt{2}\right) \\ &= \frac{3}{\sqrt{\pi}} - \frac{6 \sqrt{2}}{\pi^{3/2}}\,  \arctan (\sqrt{2}). \end{align}$$
EDIT:
Again integrating by parts twice, we get
$$\int_{0}^{\infty} \text{erfc}^{4}(x) \ dx = \frac{4}{\sqrt{\pi}} + \frac{24}{\pi} \int_{0}^{\infty} \text{erfc}^{2}(x) e^{-2x^{2}} \, dx.$$
In general, for positive parameters $a,b$, and $p$,
$$ \begin{align} I(a,b,p) &= \int_{0}^{\infty} \text{erfc}(ax) \,  \text{erfc}(bx) e^{-px^{2}} \, dx \\ &= \frac{4}{\pi} \int_{0}^{\infty} \int_{a}^{\infty} \int_{b}^{\infty} x^{2} e^{-x^{2}y^{2}} e^{-x^{2} z^{2}} e^{-px^{2}} \, dy \, dz \, dx \\ &= \frac{4}{\pi} \int_{a}^{\infty} \int_{b}^{\infty} \int_{0}^{\infty} x^{2} e^{-(p+y^{2}+z^{2})x^{2}} \, dx \, dy \, dz \\ &= \frac{1}{\sqrt{\pi}} \int_{a}^{\infty} \int_{b}^{\infty} \frac{1}{(p+y^{2}+z^{2})^{3/2}} \, dy \, dz. \end{align}$$
Let $y = \sqrt{p+z^{2}} \tan \theta$.
$$ \begin{align} I(a,b,p) &= \frac{1}{\sqrt{\pi}} \int_{a}^{\infty} \int_{\arctan(b / \sqrt{p+z^{2}})}^{\pi /2} \frac{\cos \theta}{p+z^{2}} \, d \theta \, dz \\ &= \frac{1}{\sqrt{\pi}} \int_{a}^{\infty} \frac{1}{p+z^{2}} \left(1- \frac{b}{\sqrt{p+b^{2}+z^{2}}} \right) \, dz \\ &= \frac{1}{\sqrt{\pi p}} \arctan \left(\frac{\sqrt{p}}{a} \right) - \frac{b}{\sqrt{\pi}} \int_{a}^{\infty} \frac{1}{(p+z^{2})\sqrt{p+b^{2}+z^{2}}} \, dz \end{align}$$
Now let $ \displaystyle t = \frac{1}{z}$.
$$ I(a,b,p) =\frac{1}{\sqrt{\pi p}} \arctan \left(\frac{\sqrt{p}}{a} \right) - \frac{b}{\sqrt{\pi}} \int_{0}^{1/a} \frac{t}{(1+pt^{2})\sqrt{1+b^{2}t^{2}+pt^{2}}} \, dt$$ 
Finally let $u^{2} = 1+b^{2}t^{2} + p t^{2}$.
$$ \begin{align} I(a,b,p) &= \frac{1}{\sqrt{\pi p}} \arctan \left(\frac{\sqrt{p}}{a} \right) - \frac{b}{\sqrt{\pi}} \int_{1}^{\sqrt{p+a^{2}+b^{2}}/a} \frac{1}{b^{2}+pu^{2}} \ du \\ &=\frac{1}{\sqrt{\pi p}} \arctan \left(\frac{\sqrt{p}}{a} \right) - \frac{1}{\sqrt{\pi p}}\int_{\sqrt{p} / b}^{\sqrt{p(p+a^{2}+b^{2)}}/(ab)} \frac{1}{1+w^{2}} \ dw \\ &= \frac{1}{\sqrt{\pi p}} \left[\arctan \left(\frac{\sqrt{p}}{a} \right) +\arctan \left(\frac{\sqrt{p}}{b} \right)  - \arctan \left( \frac{\sqrt{p(p+a^{2}+b^{2})}}{ab}\right) \right] \end{align}$$
Therefore,
$$ I(1,1,2) = \int_{0}^{\infty} \text{erfc}^{2}(x) e^{-2x^{2}} \ dx = \frac{1}{\sqrt{2 \pi}} \left( 2 \arctan (\sqrt{2})- \arctan (2\sqrt{2})  \right) ,$$
and
$$ \begin{align}  \int_{0}^{\infty} \text{erfc}^{4}(x) \ dx &= \frac{4}{\sqrt{\pi}} - \frac{12 \sqrt{2}}{\pi^{3/2}} \Big(2 \arctan (\sqrt{2})- \arctan (2 \sqrt{2}) \Big) \\ &= \frac{4}{\sqrt{\pi}} - \frac{24 \sqrt{2}}{\pi^{3/2}} \arctan \left( \frac{1}{2\sqrt{2}} \right). \end{align}$$
A: $I_5$ can be reduced to an integral of elementary functions. Namely,
\begin{eqnarray}
I_5&=&\int_0^{\infty}\operatorname{erfc}^5x\,dx =\\
&=&\frac{5}{\pi^{1/2}}-\frac{240\sqrt{2}\arctan{\sqrt{2}}}{\pi^{5/2}}\left(\arctan2- \frac{\pi}{4}\right)+\frac{480}{\pi^{5/2}}\int_{\operatorname{arcsinh 1}}^{\infty}\frac{\arctan(\cosh x)}{3\cosh^2x-1}dx. 
\end{eqnarray}
Making the change of variables $z=e^{-x}$ in the remaining integral, one can rewrite it as
\begin{eqnarray}
\int_{\operatorname{arcsinh 1}}^{\infty}\frac{\arctan(\cosh x)}{3\cosh^2x-1}dx=\\
=\frac{1}{\sqrt3}\Im 
\int_{1+\sqrt2}^{\infty}\left(\frac{\log\left(1+i\frac{z+z^{-1}}{2}\right)}{z^2-\frac{2}{\sqrt3}z+1}-\frac{\log\left(1+i\frac{z+z^{-1}}{2}\right)}{z^2+\frac{2}{\sqrt3}z+1}\right)dz.
\end{eqnarray}
It is rather obvious that this can be written in terms of dilogarithm values (as I suggested in the previous version of my post). In particular, Mathematica can evaluate the integrals but  I am too lazy to retype its long answer - see comment of Vladimir Reshetnikov below.

In order to obtain the first formula:


*

*Integrate twice by parts to express $I_5$ in terms of $Q=\int_0^{\infty}e^{-2x^2}\operatorname{erfc}^3x\,dx$.

*Then consider one-parameter deformation $Q(\alpha)=\int_0^{\infty}e^{-2x^2}\operatorname{erfc}^3(\alpha x)\,dx$.

*The derivative $Q'(\alpha)$ is easily computable. Integrating it back and using that $Q(\infty)=0$, we get 
\begin{eqnarray}
I_5=\frac{5}{\pi^{1/2}}-\frac{240\sqrt{2}\arctan{\sqrt{2}}}{\pi^{5/2}}\left(\arctan2- \frac{\pi}{4}\right)+\frac{480}{\pi^{5/2}}\int_1^{\infty}\frac{\arctan\sqrt{1+\alpha^2}}{\left(3\alpha^2+2\right)\sqrt{1+\alpha^2}}d\alpha. 
\end{eqnarray}

*The remaining integral gives the above expression after the change of variables $\alpha =\sinh x$.
