How to evaluate $\int_0^\infty\operatorname{erfc}^n x\ \mathrm dx$?

Question

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?

Answer 1

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.

Answer 2

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

Answer 3

Let us assume the $n \ge 2$. Then integrating by parts we have:

\begin{eqnarray} &&I_n := \int\limits_0^\infty [erfc(x)]^n dx = \\ &&\left.x \cdot [erfc(x)]^n\right|_0^\infty - \int\limits_0^\infty x \cdot n \cdot [erfc(x)]^{n-1} (-2 \sqrt{2}) \cdot \phi(\sqrt{2} x) dx =\\ &&-\sqrt{2} n \int\limits_0^\infty (\phi(\sqrt{2} x))^{'} erfc(x)^{n-1} dx = \\ &&\frac{n}{\sqrt{\pi}} - \sqrt{\frac{2}{\pi}}\cdot n(n-1) 2^{n-2} \int\limits_{-\infty}^0 \phi(x) \cdot [\Phi(\frac{x}{\sqrt{2}})]^{n-2} dx \end{eqnarray}

Now we take $n=6$ and we use use An integral involving a Gaussian and a power of a normal cumulative distribution function . this gives:

\begin{eqnarray} &&I_6 = \frac{\left(6-15 \sqrt{2}\right) \pi +120 \sqrt{2} \tan ^{-1}\left(\frac{1}{\sqrt{2}}\right)}{\pi ^{3/2}} -\\ && \frac{360 \sqrt{2}}{\pi ^{5/2}} \cdot \sum\limits_{\xi_1={-1,1}} \sum\limits_{\xi_2={-1,1}} \sum\limits_{\eta_1={-1,1}} \sum\limits_{\eta_2={-1,1}} \xi_1 \xi_2 \left({\mathfrak F}^{1,\exp(\imath \pi/4)}_{\xi_2 \sqrt{2}+\eta_2 \sqrt{3},\imath(\xi_1+\eta_1 \sqrt{2})} - {\mathfrak F}^{\exp(\imath \pi/4),\exp(\imath \pi/3)}_{\xi_2/\sqrt{2} + \eta_2 \sqrt{3/2},\xi_1\eta_1 \imath \sqrt{2+\eta_1 \sqrt{3}}} \right) \\ &&= \frac{\left(6-45 \sqrt{2}\right) \pi +120 \sqrt{2} \left(\cot ^{-1}\left(\sqrt{2}\right)+3 i \log \left(\frac{2}{4-\sqrt{2}+\sqrt{6}}\right)\right)}{\pi ^{3/2}} - \\ &&\frac{360 \sqrt{2}}{\pi ^{5/2}} \cdot \sum\limits_{\xi_1={-1,1}} \sum\limits_{\xi_2={-1,1}} \sum\limits_{\eta_1={-1,1}} \sum\limits_{\eta_2={-1,1}} \xi_1 \xi_2 \\ &&\left( \left.F[z,\xi_2 \sqrt{2}+\eta_2 \sqrt{3},\imath(\xi_1+\eta_1 \sqrt{2})]\right|_{z=1}^{z=\exp(\imath \pi/4)} - \left.F[z,\xi_2/\sqrt{2} + \eta_2 \sqrt{3/2},\xi_1\eta_1 \imath \sqrt{2+\eta_1 \sqrt{3}}]\right|_{z=\exp(\imath \pi/4)}^{z=\exp(\imath \pi/3)} \right) \end{eqnarray} Here the function ${\mathfrak F}$ is expressed through the di-logarithm and it is defined in my previous answer above.

phi[x_] := Exp[-x^2/2]/Sqrt[2 Pi]; n = 6; c = N[1/Sqrt[2], 50];
{phi2} = {ArcSin[ Sqrt[1 + 4 c^2]/Sqrt[2 + 4 c^2]]};
NIntegrate[Erfc[x]^n, {x, 0, Infinity}, WorkingPrecision -> 30]
2 Sqrt[2] n NIntegrate[
  x phi[Sqrt[2] x] Erfc[x]^(n - 1), {x, 0, Infinity}, 
  WorkingPrecision -> 30]
-Sqrt[2] n NIntegrate[
  D[phi[Sqrt[2] x], x] Erfc[x]^(n - 1), {x, 0, Infinity}, 
  WorkingPrecision -> 30]
(-(- n)/Sqrt[Pi] + (-Sqrt[2] n) (n - 1)/Sqrt[Pi] 2^(n - 2) NIntegrate[
    phi[x] (CDF[NormalDistribution[], x/Sqrt[2]])^(n - 
        2), {x, -Infinity, 0}, WorkingPrecision -> 30])
(((6 - 15 Sqrt[2]) \[Pi] + 120 Sqrt[2] ArcTan[1/Sqrt[2]])/\[Pi]^(3/2) +
  -((360 Sqrt[2])/\[Pi]^(5/2)) Sum[
    xi1 xi2 (FF[1, Exp[I Pi/4], xi2 Sqrt[2] + eta2 Sqrt[3], 
        I (xi1 + eta1 Sqrt[2])] - 
       FF[Exp[I Pi/4], Exp[I phi2], xi2/Sqrt[2] + eta2 Sqrt[3/2], 
        xi1 eta1 I Sqrt[2 + eta1 Sqrt[3]]]) , {xi2, -1, 1, 
     2}, {eta2, -1, 1, 2}, {xi1, -1, 1, 2}, {eta1, -1, 1, 2}])
N[(((6 - 45 Sqrt[2]) \[Pi] + 
    120 Sqrt[
     2] (ArcCot[Sqrt[2]] + 
       3 I Log[2/(4 - Sqrt[2] + Sqrt[6])]))/\[Pi]^(3/2) +
   -((360 Sqrt[2])/\[Pi]^(5/2)) Sum[
     xi1 xi2 (F[Exp[I Pi/4], xi2 Sqrt[2] + eta2 Sqrt[3], 
         I (xi1 + eta1 Sqrt[2])] - 
        F[1, xi2 Sqrt[2] + eta2 Sqrt[3], I (xi1 + eta1 Sqrt[2])] - 
        F[Exp[I Pi/3], xi2/Sqrt[2] + eta2 Sqrt[3/2], 
         xi1 eta1 I Sqrt[2 + eta1 Sqrt[3]]] + 
        F[Exp[I Pi/4], xi2/Sqrt[2] + eta2 Sqrt[3/2], 
         xi1 eta1 I Sqrt[2 + eta1 Sqrt[3]]]) , {xi2, -1, 1, 
      2}, {eta2, -1, 1, 2}, {xi1, -1, 1, 2}, {eta1, -1, 1, 2}]), 50]