Definite integrals over Gaussians and multiple error functions

Question

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.

Answer 1

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

Answer 2

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