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$

Question

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?

Answer 1

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.

Answer 2

Let us denote : \begin{equation} {\mathcal I}(c) := \int\limits_0^\infty \operatorname{erf}(c x) \cdot [\operatorname{erf}( x)]^2 e^{-x^2} dx \end{equation} By differentiating with respect to the parameter $c$ we have: \begin{equation} \frac{ d }{d c} {\mathcal I}(c) = \frac{2^2}{\pi^{3/2}} \frac{1}{1+c^2} \cdot \frac{1}{\sqrt{2+c^2}} \cdot \arctan\left( \frac{1}{\sqrt{2+c^2}}\right) \end{equation} therefore the only thing we need to do is to integrate the right hand side. I have calculated a more generic integral that involves this one as a special case in A generalized Ahmed's integral . Here I only state the result: \begin{eqnarray} &&{\mathcal I}(c) = \frac{4}{\pi^{3/2}} \left(\right.\\ && \arctan( \frac{c}{\sqrt{2+c^2}}) \arctan( \frac{1}{\sqrt{2+c^2}})+\\ && \frac{\imath}{2} \left.\left[ {\mathcal F}^{(\alpha_-,+e^{-\imath \phi})}(t)+ {\mathcal F}^{(\alpha_-,-e^{-\imath \phi})}(t)- {\mathcal F}^{(\alpha_-,-e^{+\imath \phi})}(t)- {\mathcal F}^{(\alpha_-,+e^{+\imath \phi})}(t) \right]\right|_0^B-\\ && \frac{\imath}{2} \left.\left[ {\mathcal F}^{(\alpha_+,+e^{-\imath \phi})}(t)+ {\mathcal F}^{(\alpha_+,-e^{-\imath \phi})}(t)- {\mathcal F}^{(\alpha_+,-e^{+\imath \phi})}(t)- {\mathcal F}^{(\alpha_+,+e^{+\imath \phi})}(t) \right]\right|_0^B \left.\right) \end{eqnarray} where $\alpha_- = \sqrt{2}-1$, $\alpha_+:=\sqrt{2}+1$, $\phi:= \arccos(1/\sqrt{3})$, $B:=(-\sqrt{2}+\sqrt{2+c^2})/c$ and \begin{eqnarray} &&{\mathcal F}^{(a,b)}(t):=\int \arctan(\frac{t}{a}) \frac{1}{t-b} dt = \log(t-b) \arctan(\frac{t}{a})\\ &&-\frac{1}{2 \imath} \left( \log(t-b) \left[ \log(\frac{t-\imath a}{b-\imath a}) - \log(\frac{t+\imath a}{b+\imath a})\right] + Li_2(\frac{b-t}{b-\imath a}) - Li_2(\frac{b-t}{b+\imath a})\right) \end{eqnarray}

Update:

Note that the anti-derivative ${\mathcal F}^{(a,b)}(t)$ may have a jump. This will happen if and only if either the quantity $(t+\imath a)/(b+\imath a)$ crosses the negative real axis or the quantity $(t-\imath a)/(b-\imath a)$ crosses the negative real axis both for some $t\in(0,B)$. This has an effect that the argument of the logarithm jumps by $2\pi$. In order to take this into account we have to exclude from the integration region a small vicinity of the singularity in question. In other words the correct formula reads:

\begin{eqnarray} &&{\mathcal I}(c) = \frac{4}{\pi^{3/2}} \left(\right.\\ && \arctan( \frac{c}{\sqrt{2+c^2}}) \arctan( \frac{1}{\sqrt{2+c^2}})+\\ && \frac{\imath}{2}\left[ {\bar {\mathcal F}}^{(\alpha_-,+e^{-\imath \phi})}(0,B)+ {\bar {\mathcal F}}^{(\alpha_-,-e^{-\imath \phi})}(0,B)- {\bar {\mathcal F}}^{(\alpha_-,-e^{+\imath \phi})}(0,B)- {\bar {\mathcal F}}^{(\alpha_-,+e^{+\imath \phi})}(0,B) \right]-\\ && \frac{\imath}{2} \left[ {\bar {\mathcal F}}^{(\alpha_+,+e^{-\imath \phi})}(0,B)+ {\bar {\mathcal F}}^{(\alpha_+,-e^{-\imath \phi})}(0,B)- {\bar {\mathcal F}}^{(\alpha_+,-e^{+\imath \phi})}(0,B)- {\bar {\mathcal F}}^{(\alpha_+,+e^{+\imath \phi})}(0,B) \right] \left.\right) \end{eqnarray}

where

\begin{eqnarray} {\bar {\mathcal F}}^{a,b}(0,B) &:=& {\mathcal F}^{(a,b)}(B)-{\mathcal F}^{(a,b)}(A) +\\ && 1_{t^{(*)}_+ \in (0,1)} \left( -{\mathcal F}^{(a,b)}(B(t^{(*)}_+ +\epsilon))+{\mathcal F}^{(a,b)}(B(t^{(*)}_+ -\epsilon))\right)+\\ && 1_{t^{(*)}_- \in (0,1)} \left( -{\mathcal F}^{(a,b)}(B(t^{(*)}_- +\epsilon))+{\mathcal F}^{(a,b)}(B(t^{(*)}_- -\epsilon))\right) \end{eqnarray}

where

\begin{eqnarray} t^{(*)}_\pm:= \frac{Im[\mp \imath a(\bar{b} \mp \imath a)]}{B Im[\bar{b} \mp \imath a]} \end{eqnarray}

See Mathematica code below for testing:

F[t_, a_, b_] := 
  Log[t - b] ArcTan[t/a] - 
   1/(2 I) (Log[
        t - b] (Log[(t - I a)/(b - I a)] - 
         Log[(t + I a)/(b + I a)]) + PolyLog[2, (b - t)/(b - I a)] - 
      PolyLog[2, (b - t)/(b + I a)]);
FF[A_, B_, a_, b_] := Module[{res, rsp, rsm, tsp, tsm, eps = 10^(-9)},
   res = F[B, a, b] - F[A, a, b];

   tsp = -(Im[I a (Conjugate[b] - I a)]/(B Im[Conjugate[b] - I a]));
   tsm = +(Im[I a (Conjugate[b] + I a)]/(B Im[Conjugate[b] + I a]));
   (*
   If[0\[LessEqual] tsp\[LessEqual]1,Print["Jump +!!"]];
   If[0\[LessEqual] tsm\[LessEqual]1,Print["Jump -!!"]];
   *)

   rsp = If[
     0 <= tsp <= 1, -F[A + (tsp + eps) (B - A), a, b] + 
      F[A + (tsp - eps) (B - A), a, b], 0];
   rsm = If[
     0 <= tsm <= 1, -F[A + (tsm + eps) (B - A), a, b] + 
      F[A + (tsm - eps) (B - A), a, b], 0];

   res + rsp + rsm

   ];

For[count = 1, count <= 100, count++,
  c = RandomReal[{-10, 10}, WorkingPrecision -> 50];
  x1 = NIntegrate[Erf[c x] Erf[x]^2 Exp[-x^2], {x, 0, Infinity}, 
    WorkingPrecision -> 30];
  4/Pi^(3/2)
    NIntegrate[
    1/(1 + xi^2) 1/Sqrt[2 + xi^2] ArcTan[1/Sqrt[2 + xi^2]], {xi, 0, 
     c}];
  A1 = 1; A2 = 1; A3 = c;
  phi = ArcCos[1/Sqrt[3]]; B = (-Sqrt[2] + Sqrt[2 + c^2])/c;
  4/Pi^(3/
    2) ((ArcTan[c/Sqrt[2 + c^2]]) ArcTan[1 /Sqrt[2 + c^2]] + 
     4  Sqrt[2]
       NIntegrate[(ArcTan[t/(Sqrt[2] - 1)] - 
          ArcTan[t/(Sqrt[2] + 1)])  
        t/((1 - t^2)^2 + (2 ) (1 + t^2)^2), {t, 
        0, (-Sqrt[2] + Sqrt[2 + c^2])/c}, WorkingPrecision -> 30]);
  4/Pi^(3/
    2) ((ArcTan[c/Sqrt[2 + c^2]]) ArcTan[1 /Sqrt[2 + c^2]] + 
     I/2 NIntegrate[(ArcTan[t/(Sqrt[2] - 1)] - 
          ArcTan[t/(Sqrt[2] + 1)]) (1/(t - E^(-I phi)) - 
          1/(t - E^(I phi)) + 1/(t + E^(-I phi)) - 
          1/(t + E^(I phi))), {t, 0, (-Sqrt[2] + Sqrt[2 + c^2])/c}, 
       WorkingPrecision -> 30]);



  x2 = 4/Pi^(3/2) ((ArcTan[c/Sqrt[2 + c^2]]) ArcTan[1/Sqrt[2 + c^2]] +
       I/2 (FF[0, B, (Sqrt[2] - 1), 1/Sqrt[3] - I Sqrt[2/3]] + 
         FF[0, B, (Sqrt[2] - 1), -(1/Sqrt[3]) + I Sqrt[2/3]] - 
         FF[0, B, (Sqrt[2] - 1), -(1/Sqrt[3]) - I Sqrt[2/3]] - 
         FF[0, B, (Sqrt[2] - 1), 1/Sqrt[3] + I Sqrt[2/3]]) - 
      I/2 (FF[0, B, (Sqrt[2] + 1), 1/Sqrt[3] - I Sqrt[2/3]] + 
         FF[0, B, (Sqrt[2] + 1), -(1/Sqrt[3]) + I Sqrt[2/3]] - 
         FF[0, B, (Sqrt[2] + 1), -(1/Sqrt[3]) - I Sqrt[2/3]] - 
         FF[0, B, (Sqrt[2] + 1), 1/Sqrt[3] + I Sqrt[2/3]]));
  If[Abs[x2/x1 - 1] > 10^(-3), 
   Print["results do not match..", {c, {x1, x2}}]; Break[]];
  If[Mod[count, 10] == 0, PrintTemporary[count]];
  ];

Answer 3

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.