An integral involving error functions and a Gaussian

Let $d\ge 1$ be an integer and let $\vec{A}:=\left\{ A_i \right\}_{i=1}^d$ be real numbers. We consider a following integral: \begin{equation} {\mathfrak I}^{(d)}(\vec{A}):=\int\limits_0^\infty e^{-u^2}\left[ \prod_{i=1}^d \operatorname{erf}(A_i u) \right] du \end{equation} By expanding the error functions in Taylor series and then integrating term by term we found the answer for $d=1$ and $d=2$. We have: \begin{eqnarray} \sqrt{\pi} {\mathfrak I}^{(d)}(\vec{A}) = \begin{cases} \arctan(A_1) & \text{if $d=1$}\\[4pt] \arctan\left(\frac{A_1 A_2}{\sqrt{1+A_1^2+A_2^2}}\right) & \text{if $d=2$} \end{cases} \end{eqnarray} Now the question is how do we derive the result for arbitrary values of $d$?

The most elegant approach to answer such questions is to establish certain recurrence relations for the quantity in question. However this task is not always straightforward to do and therefore one needs to generalize the right hand side in a clever way. Therefore we consider a more generic integral: \begin{equation} {\mathfrak I}^{(d)}_p(\vec{A}) := \int\limits_0^\infty u^p e^{-u^2} \cdot \prod\limits_{j=1}^d \operatorname{erf}(A_j u)\cdot du \end{equation} where $d\ge 0$ and $p\ge 0$ are integers. We also have ${\mathfrak I}^{(0)}_p = 1/2 ((p-1)/2)!$.

Now by differentiating the quantity above with respect to the last parameter we get a following recurrence relation: \begin{eqnarray} \frac{\partial }{\partial A_d} {\mathfrak I}^{(d)}_p(\vec{A}) = \frac{2}{\sqrt{\pi}} \cdot \frac{1}{(\sqrt{1+A_d^2})^{p+2}} \cdot {\mathfrak I}^{(d-1)}_{p+1}(\frac{\left(A_j\right)_{j=1}^{d-1}}{\sqrt{1+A_d^2}}) \end{eqnarray}

Now let us start with $d=1$. We have: \begin{eqnarray} &&{\mathfrak I}^{(1)}_p(A_1) = \frac{1}{\sqrt{\pi}} \left( \frac{p}{2}\right)! \int\limits_0^{\arctan(A_1)} \cos(\theta)^p d\theta\\ &&=\left\{ \frac{\arctan(A_1)}{\sqrt{\pi }}, \frac{A_1}{2 \sqrt{A_1^2+1}}, \frac{\left(A_1^2+1\right) \arctan(A_1)+A_1}{2 \sqrt{\pi } \left(A_1^2+1\right)}, \frac{A_1 \left(2 A_1^2+3\right)}{4 \left(A_1^2+1\right)^{3/2}},\cdots\right\} \end{eqnarray} Now we use the results above in order to derive the quantities in question for $d=2$. We have: \begin{eqnarray} &&{\mathfrak I}^{(2)}_p(\vec{A}) = \frac{2}{\sqrt{\pi}}\cdot \int\limits_0^{A_2} \frac{1}{(\sqrt{1+\xi^2})^{p+2}} \cdot {\mathfrak I}^{(1)}_{p+1}(\frac{A_1}{\sqrt{1+\xi^2}}) d\xi\\ &&=\left\{\right.\\ && \frac{1}{\sqrt{\pi}} \arctan\left(\frac{A_1 A_2}{\sqrt{1+A_1^2+A_2^2}}\right), \frac{1}{\pi}\left( \frac{A_1}{\sqrt{1+A_1^2}} \arctan(\frac{A_2}{\sqrt{1+A_1^2}}) + \frac{A_2}{\sqrt{1+A_2^2}} \arctan(\frac{A_1}{\sqrt{1+A_2^2}})\right), \frac{1}{2 \sqrt{\pi }} \left(\frac{A_1 A_2 \left(A_1^2+A_2^2+2\right) }{\left(A_1^2+1\right) \left(A_2^2+1\right) \sqrt{A_1^2+A_2^2+1}}+\arctan\left(\frac{A_1 A_2}{\sqrt{A_1^2+A_2^2+1}}\right)\right) ,\cdots\\ &&\left. \right\} \end{eqnarray} Note that the integrals which we encountered so far are doable because they involve either a product of a square root and a rational function or an arc tangent and a rational function.There are well known techniques for handling such integrals. We will provide additional results for bigger values of $d$ and $p$ later on.

Update: Now we provide the results for $d=3$. Firstly let us define: \begin{eqnarray} {\mathfrak F}^{(A,B)}_{a,b} &:=& \int\limits_A^B \frac{\log(z+a)}{z+b} dz\\ &=& F[B,a,b] - F[A,a,b] + 1_{t^* \in (0,1)} \left( -F[A+(t^*+\epsilon)(B-A),a,b] + F[A+(t^*-\epsilon)(B-A),a,b] \right) \end{eqnarray} where \begin{eqnarray} t^*:=-\frac{Im[(A+b)(b^*-a^*)]}{Im[(B-A)(b^*-a^*)]} \end{eqnarray} and \begin{equation} F[z,a,b] := \log(z+a) \log\left( \frac{z+b}{b-a}\right) + Li_2\left( \frac{z+a}{a-b}\right) \end{equation} Then the results read: \begin{eqnarray} &&{\mathfrak I}^{(3)}_p(\vec{A}) = \frac{2}{\sqrt{\pi}}\cdot \int\limits_0^{A_3} \frac{1}{(\sqrt{1+\xi^2})^{p+2}} \cdot {\mathfrak I}^{(2)}_{p+1}(\frac{(A_1,A_2)}{\sqrt{1+\xi^2}}) d\xi=\\ &&\left\{\right.\\ &&\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\! -\frac{1}{2\pi^{3/2}} \sum\limits_{\xi=1}^4 \sum\limits_{\eta=1}^4 (-1)^{\left\lfloor \frac{\eta -1}{2}\right\rfloor +\left\lfloor \frac{\xi -1}{2}\right\rfloor }\cdot \left({\mathfrak F}^{(1,\frac{\sqrt{1+A_1^2+A_3^2} - |A_3|}{\sqrt{1+A_1^2}})}_{\frac{i \left((-1)^{\xi } \sqrt{A_1^2+A_2^2+1}+A_2 (-1)^{\left\lfloor \frac{\xi -1}{2}\right\rfloor }\right)}{\sqrt{A_1^2+1}},\frac{A_1 (-1)^{\eta }+i (-1)^{\left\lceil \frac{\eta -1}{2}\right\rceil }}{\sqrt{A_1^2+1}}} + % {\mathfrak F}^{(1,\frac{\sqrt{1+A_2^2+A_3^2} - |A_3|}{\sqrt{1+A_2^2}})}_{\frac{i \left((-1)^{\xi } \sqrt{A_2^2+A_1^2+1}+A_1 (-1)^{\left\lfloor \frac{\xi -1}{2}\right\rfloor }\right)}{\sqrt{A_2^2+1}},\frac{A_2 (-1)^{\eta }+i (-1)^{\left\lceil \frac{\eta -1}{2}\right\rceil }}{\sqrt{A_2^2+1}}} \right),\\ &&\frac{1}{\pi} \sum\limits_{j=1}^3 \frac{A_j}{\sqrt{1+A_j^2}}\cdot \arctan\left( \frac{\prod\limits_{l=1,l\ne j}^3 A_l}{\sqrt{1+A_j^2} \sqrt{1+A_1^2+A_2^2+A_3^2}}\right),\\ &&\cdots \\ &&\left. \right\} \end{eqnarray} Now let us take $d=4$. Let us define: \begin{eqnarray} S&:=&\sum\limits_{l=1}^3 A_l^2 \end{eqnarray} and $l\%2 = 1_{l\in {\mathbb N} \setminus 2{\mathbb N}} + 2 \cdot 1_{l\in 2{\mathbb N}}$. Then we define: \begin{eqnarray} \bar{r}_{j,l_1} &:=& \imath \frac{(-1)^{\lfloor \frac{l_1-1}{2}\rfloor} \sqrt{A_j^2(1+S)} + (-1)^{l_1} \sqrt{S(1+A_j^2)}}{\sqrt{S-A_j^2}} \quad \mbox{for $j=1,\cdots,3$ and $l_1=1,\cdots,4$}\\ r_{j,l}&:=& \frac{(-1)^{1+\lfloor \frac{l-1}{2}\rfloor}A_{m_{j,l\%2}} \sqrt{1+S)} + (-1)^{1+\lfloor \frac{l-1}{4}\rfloor} \imath \sqrt{(1+A_j^2)(S-A_j^2-A_{m_{j,l\%2}}^2)}}{\sqrt{(1+A_j^2+A_{m_{j,l\%2}}^2)(S-A_j^2)}} \quad \mbox{for $j=1,\cdots,3$ and $l=1,\cdots,8$}\\ \end{eqnarray}

And now we have: \begin{eqnarray} &&{\mathfrak I}^{(4)}_p(\vec{A}) = \frac{2}{\sqrt{\pi}}\cdot \int\limits_0^{A_4} \frac{1}{(\sqrt{1+\xi^2})^{p+2}} \cdot {\mathfrak I}^{(3)}_{p+1}(\frac{(A_1,A_2,A_3)}{\sqrt{1+\xi^2}}) d\xi=\\ &&\left\{\right.\\ % &&\frac{2}{\pi^{3/2}} \sum\limits_{j=1}^3 \arctan\left( \frac{A_j A_4}{\sqrt{1+A_j^2+A_4^2}}\right) \arctan\left( \frac{\prod\limits_{l=1,l\ne j}^3 A_l}{\sqrt{1+A_j^2+A_4^2} \sqrt{1+S+A_4^2}}\right)+\\ &&\frac{1}{2 \pi^{3/2}} \sum\limits_{l=1}^8 \sum\limits_{l_1=1}^4 \sum\limits_{j=1}^3 (-1)^{\left\lfloor \frac{l-1}{4}\right\rfloor +\left\lfloor \frac{l-1}{2}\right\rfloor +\left\lfloor \frac{l_1-1}{2}\right\rfloor} {\mathfrak F}^{(0,\frac{\sqrt{(1+S)(1+A_j^2+A_4^2)} - \sqrt{(1+A_j^2)(1+S+A_4^2)}}{A_4 \sqrt{S-A_j^2}})}_{-\bar{r}_{j,l_1},-r_{j,l}},\\ &&\cdots\\ &&\left.\right\} \end{eqnarray}