# 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?

• Perhaps include how you did the first two, someone may have luck generalising. – Bennett Gardiner Mar 2 '14 at 20:37
• In Prudnikov-Brychkov-Marychev, vol. 2 there is a formula for $$\int_0^{\infty} \prod_{k=1}^4 \operatorname{erfc}(c_k x)\,dx$$ in terms of elementary functions of parameters $c_1\ldots c_4$. – Start wearing purple Mar 14 '14 at 20:09

Here is my trial, which is partially successful but still not fully answering to your question.

Using some coordinate change, I derived that

$$I_{n} := \int_{0}^{\infty} \mathrm{erfc}^{n}(x) \, dx = \frac{1}{\sqrt{n}} \left( \frac{2}{\sqrt{\pi}} \right)^{n-1} \int_{T^{n-1}} \int_{0}^{\infty} s^{n-1}e^{-(|x|^{2}-1)s^{2}} \mathrm{erfc}(s) \, ds d\sigma_{x}, \tag{*}$$

where $T^{n-1} = \{ x \in [0, \infty)^{n} : x_{1} + \cdots + x_{n} = \sqrt{n} \}$ is an $(n-1)$-simplex and $d\sigma_{x}$ is the surface measure on $T^{n-1}$.

Case $n = 1$. Since $T^{0} = \{1\}$ and $d\sigma_{x} = \delta(x-1)$ is a unit mass located at $x = 1$, the equation $\text{(*)}$ gives no new information.

Case $n = 2.$ It is easy to find that

$$\int_{0}^{\infty} s e^{-(|x|^{2}-1)s^{2}} \mathrm{erfc}(s) \, ds = \frac{1}{2|x|(|x|+1)}$$

and by suitable parametrization of the line segment $T^{1}$, we get

\begin{align*} I_{2} &= \sqrt{\frac{2}{\pi}} \int_{T^{1}} \frac{1}{2|x|(|x|+1)} \, d\sigma_{x} \\ &= \sqrt{\frac{2}{\pi}} \int_{0}^{\sqrt{2}} \frac{\sqrt{2} \, dt}{2\sqrt{t^{2} + (\sqrt{2}-t)^{2}} ( \sqrt{t^{2} + (\sqrt{2}-t)^{2}} + 1)} \\ &= \sqrt{\frac{2}{\pi}} \int_{0}^{\frac{\pi}{4}} \frac{d\theta}{1+\cos\theta} = \sqrt{\frac{2}{\pi}} \tan\left(\frac{\pi}{8}\right) = \frac{2 - \sqrt{2}}{\sqrt{\pi}}. \end{align*}

We can also write

$$I_{2} = \boxed{ \displaystyle \frac{2}{\sqrt{\pi}} - \frac{2\sqrt{2}}{\pi^{3/2}} \arctan(\infty) }.$$

Case $n = 3$. We have

$$\int_{0}^{\infty} s^{2}e^{-(|x|^{2}-1)s^{2}} \mathrm{erfc}(s) \, ds = \frac{1}{2\sqrt{\pi}(|x|^{2}-1)} \left( \frac{\arctan\sqrt{|x|^{2}-1}}{\sqrt{|x|^{2}-1}} - \frac{1}{|x|^{2}} \right).$$

Now we know that $T^{2}$ is a regular triangle with side length $\sqrt{6}$. Thus by introducing polar coordinate change, we get

\begin{align*} I_{3} &= \frac{2}{\pi^{3/2}\sqrt{3}} \int_{T^{2}} \frac{1}{|x|^{2}-1} \left( \frac{\arctan\sqrt{|x|^{2}-1}}{\sqrt{|x|^{2}-1}} - \frac{1}{|x|^{2}} \right) \, d\sigma_{x} \\ &= \frac{4\sqrt{3}}{\pi^{3/2}} \int_{0}^{\frac{\pi}{3}} \int_{0}^{\frac{\sec\theta}{\sqrt{2}}} \frac{1}{r} \left( \frac{\arctan r}{r} - \frac{1}{r^{2} + 1} \right) \, drd\theta \\ &= \frac{4\sqrt{3}}{\pi^{3/2}} \int_{0}^{\frac{\pi}{3}} \left( 1- \frac{\arctan\left(\sec\theta / \sqrt{2}\right)}{\sec\theta / \sqrt{2}} \right) \, d\theta \\ &= \frac{4\sqrt{3}}{\pi^{3/2}} \int_{\frac{1}{\sqrt{2}}}^{\sqrt{2}} \left( 1 - \frac{\arctan u}{u} \right) \frac{du}{u\sqrt{2u^{2} - 1}}. \end{align*}

I haven't tried the last integral, but Mathematica suggests that

$$\int_{\frac{1}{\sqrt{2}}}^{\sqrt{2}} \left( 1 - \frac{\arctan u}{u} \right) \frac{du}{u\sqrt{2u^{2} - 1}} = \frac{\sqrt{3}\pi}{4} -\sqrt{\frac{3}{2}} \arctan\left( \sqrt{2} \right).$$

This gives

$$I_{3} = \boxed{ \displaystyle \frac{3}{\sqrt{\pi}} - \frac{6\sqrt{2}}{\pi^{3/2}} \arctan\left( \sqrt{2} \right) },$$

which can be numerically checked.

I am trying to generalize this calculation for higher dimensions, but it seems somewhat daunting. For example, if $n = 4$ we have to evaluate

$$\int_{0}^{\infty} \mathrm{erfc}^{4}(x) \, dx = \frac{1}{\pi^{3/2}} \int_{T^{3}} \frac{2|x|+1}{(|x|+1)^{2}|x|^{3}} \, d\sigma_{x}.$$

Nevertheless, using the formula

$$\int_{T^{n-1}} f(|x|) \, d\sigma_{x} = n! \int_{0}^{\sqrt{\frac{1}{n-1}}} \int_{0}^{\sqrt{\frac{n}{n-2}}s_{1}} \cdots \int_{0}^{\sqrt{\frac{3}{1}}s_{n-2}} f\left( \sqrt{1+s_{1}^{2} + \cdots + s_{n-1}^{2}} \right) \, ds_{n-1} \cdots ds_{1}$$

together with aid of Mathematica, I was able to evaluate $I_{4}$ and it was

$$I_{4} = \boxed{ \displaystyle \frac{4}{\sqrt{\pi}} - \frac{24\sqrt{2}}{\pi^{3/2}} \arctan \left( \frac{1}{\sqrt{8}} \right) }.$$

These series of observations lead us to believe that $I_{n}$ is of the form

$$I_{n} = \frac{n}{\sqrt{\pi}} - \frac{n!\sqrt{2}}{\pi^{3/2}} \arctan \alpha_{n}$$

for some reasonably nice $\alpha_{n}$ (with $\alpha_{1} = 0$, $\alpha_{2} = \infty$, $\alpha_{3} = \sqrt{2}$ and $\alpha_{4} = \frac{1}{\sqrt{8}}$), but inverse symbolic calculations show that this seems no longer the case for $n \geq 5$.

Indeed, for $n = 5$ we have

$$I_{5} = \frac{5}{\sqrt{\pi}} - \frac{80\sqrt{2}}{\pi^{3/2}} \arctan \sqrt{\frac{2}{3}} + \frac{240\sqrt{2}}{\pi^{5/2}} A\left( \sqrt{\frac{5}{3}}, \sqrt{\frac{3}{2}}, \sqrt{\frac{4}{5}} \right),$$

where $A(p, q, r)$ is the generalized Ahmed's integral. I have no idea whether it will reduce to a familiar closed form or not.

• Very nice attempt! – Bennett Gardiner Mar 4 '14 at 22:48
• I suggest +100 bounty for a closed-form solution for $n=5$. – Vladimir Reshetnikov Mar 8 '14 at 19:54
• could you please explain you first step, When you derived your first Equation. – adam Feb 27 '19 at 10:23

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}

• Incredibly simple! A fantastic calculation + 1 upvote. – Sangchul Lee Mar 4 '14 at 23:12
• Thanks. I managed to work out something similar (but not quite as simple) for $n=4$. Things seem to get a lot more complicated for $n=5$, though. – Random Variable Mar 4 '14 at 23:37
• Yeah, your approach reduces the degree of $\mathrm{erfc}$ by 2, so in $n = 5$ we will have something involving $\int_{0}^{\infty} \mathrm{erfc}^{3}(x) e^{-2x^{2}} \, dx$. It seems still daunting... – Sangchul Lee Mar 5 '14 at 0:13
• @sos440 I'm having quite a difficult time trying to evaluate $\int_{0}^{\infty} \text{erfc}(ax) \text{erfc}(bx) \text{erfc}(cx) e^{-px^{2}} \ dx$. But interestingly $\int_{0}^{\infty} x \ \text{erf}(ax) \text{erf}(bx) \text{erf}(cx) e^{-px^{2}} \ dx$ does have a nice closed-form expression. mathworld.wolfram.com/Erf.html (35) – Random Variable Mar 5 '14 at 7:07
• Wow, that is sort of unexpected! – Sangchul Lee Mar 5 '14 at 8:31

$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:

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

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

3. 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}

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

• Thanks! The result is 5/√π^5(π(π+9√2π-6√2ArcTan[(1240029+2590550√2)/1763303])+3√2(8ArcTan[3/4]ArcTan[3-√8]-16ArcTan[2]ArcTan[√2]+ArcTan[(4√2-3)(2+3√3)/23]^2-4ArcTan[√2-√3-√6]^2))+120√2/√π^5Re[-PolyLog[2,(9+6√2)/(2I√2-1)]/2+PolyLog[2,I((8+3I)+4√(4+3I))/3]/2+PolyLog[2,(3-3I)(I-√2)/(3-3√2+√(3-6I√2))]+PolyLog[2,(3-3I)(√2-I)/(3√2-3+√(3-6I√2))]+PolyLog[2,(3+3I)(1+√2)/(3+3√2+√(3-6I√2))]-PolyLog[2,(3-3I)(I-√2)/(3-3√2+√(3+6I√2))]-PolyLog[2,(3-3I)(√2-I)/(3√2-3+√(3+6I√2))]-PolyLog[2,(3+3I)(1+√2)/(3+3√2+√(3+6I√2))]-PolyLog[2,(3+3I)(1+√2)/(3+3√2-√(3+6I√2))]+PolyLog[2,(3+3I)(1+√2)/(3+3√2-√(3-6I√2))]] – Vladimir Reshetnikov Mar 16 '14 at 0:43
• Hopefully, PolyLog part can be simplified. – Vladimir Reshetnikov Mar 16 '14 at 0:44

Let us assume $$n>3$$. We have: $$\begin{eqnarray} &&\int\limits_0^\infty \mbox{erfc}(x)^n dx =\\ && -\frac{2^n}{\sqrt{2}} \int\limits_{(-\infty,0]^n} \phi(\xi_1) \cdot \dots \cdot \phi(\xi_n) \cdot \left( \xi_1 \vee \dots \vee \xi_n\right) \cdot \prod\limits_{j=1}^n d\xi_j=\\ && -\frac{2^n}{\sqrt{2}} n! \cdot \int\limits_{-\infty < \xi_1<\xi_2<\dots<\xi_{n-1}<0} \left(\prod\limits_{j=1}^{n-1} \phi(\xi_j)\right) \cdot \left[\phi(\xi_{n-1})-\phi(0) \right]\cdot \prod\limits_{j=1}^{n-1} d\xi_j=\\ && -\frac{2^n}{\sqrt{2}} n! \cdot \left[ -\frac{1}{4 \sqrt{\pi}} \frac{\Phi(0)^{n-2}}{(n-2)!}-\phi(0) \frac{\Phi(0)^{n-1}}{(n-1)!} + \right.\\ &&\left. \frac{1}{2 \sqrt{\pi}} \int\limits_{-\infty<\xi_1 < \dots < \xi_{n-2}<0 }\left( \prod\limits_{j=1}^{n-3} \phi(\xi_j)\right) \cdot \phi(\xi_{n-2}) \Phi(-\sqrt{2} \xi_{n-2}) \prod\limits_{j=1}^{n-2} d\xi_j \right]=\\ && -\frac{2^n}{\sqrt{2}} n! \cdot \left[ \frac{2^{1-n} \arctan\left(\sqrt{2}\right)}{\pi ^{3/2} (n-3)!}-\frac{2^{\frac{1}{2}-n}}{\sqrt{\pi } (n-1)!} + \right.\\ &&\left.\frac{1}{2\sqrt{\pi}} \int\limits_{-\infty<\xi_1< \dots < \xi_{n-3}<0} \left( \prod\limits_{j=1}^{n-3} \phi(\xi_j) \right) \cdot \left[ T(\xi_{n-3},-\sqrt{2})\right] \prod\limits_{j=1}^{n-3} d\xi_j \right]\\ %&&\frac{n(n+\sqrt{2}-1)}{\sqrt{2 \pi}} - \frac{2^{n-1} n!}{\sqrt{2 \pi}} %\int\limits_{-\infty}^0 \frac{[\Phi(\xi)]^{n-3}}{(n-3)!} \cdot \Phi(-\sqrt{2} %\xi) \frac{e^{-\frac{1}{2} \xi^2}}{\sqrt{2 \pi}} d\xi \end{eqnarray}$$ In the first line we expressed the integrand through the CDF and then we used its integral definition and we integrated over $$x$$. In the second line we integrated over $$\xi_n$$ by using the identity $$\phi^{'}(\xi)=-\xi \phi(\xi)$$. In the third line we extracted the term proportional to $$\phi(0)$$ and in the remaining term we integrated over $$\xi_{n-1}$$ by completing the relevant terms to square. In the fourth line we integrated over $$\xi_{n-2}$$ by using known Gaussian integrals (see here https://en.wikipedia.org/wiki/List_of_integrals_of_Gaussian_functions for example).

Now let us assume $$n=5$$. Then the remaining integral in the bottom line above is two dimensional. We have: $$\begin{eqnarray} &&rhs=\\ &&-\frac{2^n}{\sqrt{2}} n! \cdot \left[ -\frac{4+\sqrt{2}}{768 \sqrt{\pi}} + \frac{1}{2\sqrt{\pi}} \int\limits_{-\infty}^0 \Phi(-\sqrt{2} \xi_2) [\Phi(\xi_2)]^2 \frac{1}{2} \phi(\xi_2) d\xi_2 \right]=\\ &&-\frac{2^n}{\sqrt{2}} n! \cdot \left[ \frac{\left(-28-\sqrt{2}\right) \pi +72 \tan ^{-1}\left(\sqrt{2}\right)}{768 \pi ^{3/2}} + \frac{1}{2\sqrt{\pi}} \int\limits_{-\sqrt{2}}^\infty \frac{\arctan(\frac{1}{\sqrt{1+c^2}})}{4\pi^2 (1+c^2)\sqrt{2+c^2}}dc \right]=\\ &&-\frac{2^n}{\sqrt{2}} n! \cdot \left[ -\frac{1}{12 \sqrt{\pi }}-\frac{1}{384 \sqrt{2 \pi }}+\frac{5 \tan ^{-1}\left(\sqrt{2}\right)}{32 \pi ^{3/2}}+\right. \left. \frac{\imath}{8 \pi^{3/2}}\left(-2\int\limits_{\exp(\imath \pi/4)}^{\exp(\imath \pi/2)}\frac{\left(z^2-1\right) \left(\log \left(z^2+2 \sqrt{2} z-1\right)-\log \left(z^2-2 \sqrt{2} z-1\right)\right)}{z^4+6 z^2+1} dz+ \int\limits_{\exp(\imath \pi/4)}^{1}\frac{\left(z^2-1\right) \left(\log \left(z^2+2 \sqrt{2} z-1\right)-\log \left(z^2-2 \sqrt{2} z-1\right)\right)}{z^4+6 z^2+1} dz \right) \right] \end{eqnarray}$$ In the first line we have integrated over $$\xi_1$$ and then over $$\xi_2$$ by parts once. In the second line we replaced $$-\sqrt{2} \rightarrow c$$ then differentiated with respect to $$c$$ then collected to square the $$\xi_2$$-powers in the exponentials and integrated over $$\xi_2$$ and finally we integrated over $$c$$ back again subject to the boundary values at plus infinity being equal to zero. In the third line we firstly substituted for $$1/\sqrt{2+c^2}$$ then we used an appropriate trigonometric substitution and finally we substituted for the phase factor (the exponential of the angle times the imaginary unit).Now, the remaining integrals can always be reduced to di-logarithms. Unfortunately I don't have time to complete this mundane task now and will try to finish it asap.

Update: Let us define: $$\begin{eqnarray} {\mathfrak F}^{(A,B)}_{a,b} &:=& \int\limits_A^B \frac{\log(z+a)}{z+b} dz\\ &=& F[B,a,b] - F[A,a,b] + 1_{t^* \in (0,1)} \left( -F[A+(t^*+\epsilon)(B-A),a,b] + F[A+(t^*-\epsilon)(B-A),a,b] \right) \end{eqnarray}$$ where $$\begin{eqnarray} t^*:=-\frac{Im[(A+b)(b^*-a^*)]}{Im[(B-A)(b^*-a^*)]} \end{eqnarray}$$ and $$$$F[z,a,b] := \log(z+a) \log\left( \frac{z+b}{b-a}\right) + Li_2\left( \frac{z+a}{a-b}\right)$$$$ Then the result reads: $$\begin{eqnarray} &&\int\limits_0^\infty \mbox{erfc}(x)^5 dx = \\ &&-\frac{2^5}{\sqrt{2}} 5! \left(\right.\\ &&\left. -\frac{\left(\sqrt{2}-20\right) \pi +12 \tan ^{-1}\left(2 \sqrt{2}\right)}{768 \pi ^{3/2}}+\right.\\ &&\left. \frac{1}{64 \sqrt{2} \pi^{5/2}} \sum\limits_{\xi \in \{-1,1\}} \sum\limits_{\eta \in \{-1,1\}} \sum\limits_{\xi_1 \in \{-1,1\}} \sum\limits_{\eta_1 \in \{-1,1\}} % \xi \xi_1\frac{4 \eta + 2 \sqrt{2}}{\sqrt{3+2 \eta \sqrt{2}}} \left( 2{\mathfrak F}^{(e^{\imath \pi/4},e^{\imath \pi/2})}_{\xi_1 \sqrt{2}+\eta_1 \sqrt{3}, -\imath \xi \sqrt{3+2 \eta \sqrt{2}}}- {\mathfrak F}^{(e^{\imath \pi/4},1)}_{\xi_1 \sqrt{2}+\eta_1 \sqrt{3}, -\imath \xi \sqrt{3+2 \eta \sqrt{2}}} \right) \right.\\ &&\left. \right) \end{eqnarray}$$

n = 5;
I1 = NIntegrate[Erfc[x]^n, {x, 0, Infinity}, WorkingPrecision -> 30]

N[-2^n/Sqrt[
2] 5! (-(((-20 + Sqrt[2]) \[Pi] + 12 ArcTan[2 Sqrt[2]])/(
768 \[Pi]^(3/2))) +
1/(64 Sqrt[2] \[Pi]^(5/2))
Sum[((4 eta + 2 Sqrt[2]) xi)/Sqrt[3 + 2 eta Sqrt[2]]
xi1 (2 FF[E^(I Pi/4), E^(I Pi/2),
xi1 Sqrt[2] + eta1 Sqrt[3], -I xi Sqrt[
3 + 2 eta Sqrt[2]]] -
FF[E^(I Pi/4), 1,
xi1 Sqrt[2] + eta1 Sqrt[3], -I xi Sqrt[
3 + 2 eta Sqrt[2]]]), {eta, -1, 1, 2}, {xi, -1, 1,
2}, {eta1, -1, 1, 2}, {xi1, -1, 1, 2}]), 50]