# Definite integrals over Gaussians and multiple error functions

Mathematica will gladly tell me that the integral

$$I\left[y,a\right]=\int_{y}^{\infty}dx\,e^{-x^{2}}\mathrm{erf}\left(ax\right)$$

where $$\mathrm{erf}(x)=\frac{2}{\sqrt{\pi}}\int_{0}^{x}dt\,e^{-t^{2}}$$ is the error function, can be written as

$$I\left[y,a\right]=-\frac{1}{2} \sqrt{\pi } \left(4 T\left(\sqrt{2} a y,\frac{1}{a}\right)+\mathrm{erf}(y)\, \mathrm{erf}(a y)-1\right)$$

where $$T(x,a) =\frac{1}{2 \pi }\int_0^a \frac{e^{-\left(t^2+1\right) x^2/2}}{ t^2+1} \, dt$$ is Owens T-function.

How is this derived? And more importantly: Can a similar result be derived for multiple error functions like $$I\left[y;a_{1},\ldots a_{n}\right]=\int_{y}^{\infty}dx\,e^{-x^{2}}\prod_{j=1}^{n}\mathrm{erf}\left(a_{j}x\right)$$

My end goal is to compute integrals like $$G_{n}=\int_{-\infty}^{\infty}dx_{0}\,\prod_{j=1}^{n}\int_{x_{j-1}}^{\infty}dx_{j}\,e^{-\sum_{j=0}^{n}x_{j}^{2}} \prod_{j=0}^n x_j^{p_j}$$ where the $$p_j$$ are "not-too-large" non-negative integers. In these integrals, multiple error functions naturally pop up.

• I think the result follows from differentiation with respect to the integral sign w.r.t. $a$. Apr 29, 2021 at 13:54
• @projectilemotion Thank you. I now also realize that mathematica is (as usual?) not giving me the "simplest" for of the T-function. Apr 29, 2021 at 14:15
• Indeed, I obtain the result $2\sqrt{\pi}\cdot T(\sqrt{2}y,a)$, which agrees numerically with the result you wrote. Apr 29, 2021 at 14:41
• Also, I don't think a "nice" result for multiple error functions can be derived for multiple error functions. Even the relatively simple case $y=0$ has complicated results for $n\geq 2$: An integral involving error functions and a Gaussian. Apr 29, 2021 at 15:23
• @projectilemotion Indeed it get's pretty ugly very quickly. It seems however that one can keep the number or error functions down by converting them to T-functions, whether that helps in the long run remains to be seen... Apr 30, 2021 at 8:35

This is not a full answer, more of a note. Using integration by parts, notice that: $$\int_0^\infty e^{-x^2}\operatorname{erf}(ax)\,dx=\frac{\sqrt{\pi}}{2}-\int_0^\infty e^{-x^2}\operatorname{erf}\left(\frac xa\right)\,dx$$ and you can try to split your integral up into: $$\int_y^\infty=\int_0^\infty-\int_0^y$$

Addressing what others have said: $$I(y,a)=\int_y^\infty e^{-x^2}\operatorname{erf}(ax)\,dx$$ $$\frac{\partial I(y,a)}{\partial a}=\int_y^\infty e^{-x^2}\frac{\partial}{\partial a}\left[\operatorname{erf}(ax)\right]\,dx=\frac{2}{\sqrt{\pi}}\int_y^\infty xe^{-(a^2+1)x^2}dx$$ $$=\frac{1}{\sqrt{\pi}}\frac{e^{-(a^2+1)y^2}}{(a^2+1)}$$ now you need to integrate wrt $$a$$, which you will see brings in our $$T$$ function

• Thanks for noticing. The partial integration does however not really do the trick since that just transform my integral into a integral of the same form , but with a different $y$ and $a$. Apr 29, 2021 at 14:05
• true but the first part of your question was where did the $T$ function come from, and that is the answer Apr 29, 2021 at 19:01
• I acknowledge that the comment's about integration under the integral sign got me on the right track. I had been spending the entire day playing the various forms of partial integration and not getting anywhere... Apr 30, 2021 at 8:33
• This comment was indeed very helpful and allowed me to write a compute algebra script to do the integrals for me. You can read about it here: arxiv.org/abs/2202.01090 Feb 21, 2023 at 16:24
• @MikaelFremling glad to hear it! Feb 25, 2023 at 5:20

Let $$n \ge 1$$ be an integer and let $$\vec{a} \in {\mathbb R}^n$$ and $$y \in {\mathbb R}_+$$. Then the integral in question can be thought of as a vector argument Owen's T function. In other words we have:

$$\begin{eqnarray} T[h,\vec{a}] &:=& \int\limits_y^\infty \left(\prod\limits_{i=1}^n \frac{1}{2}\operatorname{erf} [ \frac{a_j x}{\sqrt{2}}]\right) \cdot \frac{e^{-\frac{x^2}{2}}}{\sqrt{2\pi}} dx \\ &=&\frac{1}{2 \pi^{\frac{n+1}{2}}} \int\limits_{\otimes_{j=1}^n [0, a_j]} \frac{\Gamma[\frac{n+1}{2}, \frac{y^2}{2} (1+\sum\limits_{j=1}^n z_\xi^2)]}{\sqrt{1 + \sum\limits_{j=1}^n z_\xi^2}^{n+1}} \cdot \prod\limits_{j=1}^n dz_\xi \tag{1} \end{eqnarray}$$

In[563]:=
n = RandomInteger[{2, 4}];
Clear[a]; z =.; x =.; j =.;
y = RandomReal[{0, 2}];
Do[ a[j] = RandomReal[{0, 1}], {j, 1, n}];
NIntegrate[
Product[1/2 Erf[a[j] x/Sqrt[2]], {j, 1, n}] Exp[-x^2/2]/
Sqrt[2 Pi], {x, y, Infinity}]

1/2 \[Pi]^(-(1/2) - n/2)
NIntegrate[ Gamma[(n + 1)/2, y^2/2 (1 + Sum[z[xi]^2, {xi, 1, n}])]/
Sqrt[1 + Sum[z[xi]^2, {xi, 1, n}]]^(n + 1),
Evaluate[Sequence @@ Table[{z[xi], 0, a[xi]}, {xi, 1, n}]]]

Out[567]= 0.000468798

Out[568]= 0.000468798


Clearly as $$n= 1$$ the result reduces to the ordinary Owen's T function as in Wikipedia.

But now, the question appears can we come up with some fast and efficient numerical algorithm for evaluating that function? Could we generalize this algorithm when $$n >1$$?

• Thanks for the input. In the meantime i wrote an algebra script to do it for me. You can read about it here: arxiv.org/abs/2202.01090 Feb 21, 2023 at 16:22