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

• have you tried differentiation under the integral sign? the resulting integrallooks doable using integration by parts...maybe integrating back w.r.t. $c$ also isn't that bad – tired Feb 4 '17 at 15:51

Let us denote : $$$${\mathcal I}(c) := \int\limits_0^\infty \operatorname{erf}(c x) \cdot [\operatorname{erf}( x)]^2 e^{-x^2} dx$$$$ By differentiating with respect to the parameter $$c$$ we have: $$$$\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)$$$$ 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]];
];

• Thanks for answer,but yours formula holds only for -2.7 < c < 2.7. – Mariusz Iwaniuk Jul 20 '17 at 20:20
• @MariuszIwaniuk: Yes you are right something weird happens above the magical values $|c| > 2.61$. I think that the formula itself is correct there is only a problem with implementing it in computer software. Set $a=\sqrt{2}-1$, $b= -1/\sqrt{3} + \imath \sqrt{2/3}$, $t=(-\sqrt{2} + \sqrt{2+c^2})/c$ and plot both the real and the imaginary part of $Li_2((b-t)/(b-\imath a)$ as a function of $c$. What you will see is that at that magical value the real part has a kink and the imaginary part a jump which is of course wrong. There should be continuous dependance on $c$. – Przemo Jul 21 '17 at 10:43
• ,Maybe You can improve yours formula.There will be a reward :) – Mariusz Iwaniuk Jul 21 '17 at 16:03
• @MariuszIwaniuk: Well the formula is correct. If you want to understand how it was derived go to my post linked above. I am just saying that the reason why the formula ceases to work for "big" values of $c$ is that computing dilogarithms is badly implemented in the computer software that you were using. Note that we are outside the radius of convergence so we need to use analytic continuation. – Przemo Jul 21 '17 at 16:44
• @ MariuszIwaniuk: I was wrong in my comment above. I flagellate myself ;-). The implementation of the di-logarithm in Mathematica is already correct and the reason why the formula did not work for all real $c$ is explained in my update above. now the new formula does work for all real values of $c$. Feel free to use the code above to verify that. Can I get a reward now please;-). – Przemo Apr 11 '19 at 17:25

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.

• Jack, have you looked at the result after integrating over just one of $a,b,c$? It looks more than ugly - Medusa-like; avert your eyes. – Mark Viola Feb 4 '17 at 16:08
• even the second integration becomes a total mess (in my eyes) ...are you sure this works out? btw. my approach which only differentiates once is also doomed because backintegration seems impossible – tired Feb 4 '17 at 16:11
• @Dr.MV: one integration is doable by hand, the following two are not. However, we may adjust the integration domain for the $b,c$ variables (making it a quarter of circle instead of a square) and getting pretty good approximations. – Jack D'Aurizio Feb 4 '17 at 16:11
• Just curious, but was the goal to obtain a good approximation? – Mark Viola Feb 4 '17 at 16:12
• @Dr.MV: that is a question only the OP may answer. I just tried to do my best. – Jack D'Aurizio Feb 4 '17 at 16:15

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.

• H_{0,3}(c) is essentially the famous Ahmed integral so at least for some special cases like $c=1$ i think we can expect a closed form solution mathworld.wolfram.com/AhmedsIntegral.html – tired Feb 5 '17 at 9:15
• @tired You're right, I had not thought about the Ahmed's integral! I remember a paper with some generalization about this integral,maybe it contains something interesting. – Marco Cantarini Feb 5 '17 at 11:36