An integral involving a Gaussian, error functions and the Owen's T function.

This question is closely related to An integral involving a Gaussian and an Owen's T function. and An integral involving error functions and a Gaussian .

Let $$\nu_1 \ge 1$$ and $$\nu_2 \ge 1$$ be integers . Now let $$\vec{a}:=\left\{a_i\right\}_{i=1}^{\nu_1}$$ and $$\vec{b}:=\left\{b_i\right\}_{i=1}^{\nu_2}$$ and $$\vec{c}:=\left\{ c_i\right\}_{i=1}^{\nu_2}$$ be real numbers, let $$\rho(\xi):=\exp(-1/2 \xi^2)/\sqrt{2\pi}$$ and let $$T(\cdot,\cdot)$$ be the Owen's T function. Consider a following integral: $$$${\mathfrak J}^{(\nu_1,\nu_2)}(\vec{a},\vec{b},\vec{c}) := \int\limits_0^\infty \rho(\xi) \left[\prod\limits_{i=1}^{\nu_1} erf(a_i \xi)\right] \cdot \left[ \prod\limits_{i=1}^{\nu_2} T(b_i \xi, c_i)\right] d\xi \quad (i)$$$$

Now we managed to find a closed form solution for the integral above in case $$\nu_1=\nu_2=1$$. 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 $$$$F[z,a,b] := \log(z+a) \log\left( \frac{z+b}{b-a}\right) + Li_2\left( \frac{z+a}{a-b}\right)$$$$ for $$a$$,$$b$$,$$A$$,$$B$$ being complex.

Then we have: $$\begin{eqnarray} &&{\mathfrak J}^{(1,1)}(a,b,c)= \frac{1}{\pi^2}\cdot \left(\right.\\ &&\left. \frac{1}{2} \arctan(\sqrt{2} a) \arctan(c) - \frac{1}{8} \sum\limits_{i=1}^4 \sum\limits_{j=1}^4 (-1)^{j-1+\lfloor \frac{i-1}{2} \rfloor } % {\mathfrak F}^{(1,\frac{\sqrt{1+2 a^2+b^2} - \sqrt{2} a}{\sqrt{1+b^2}})}_{\frac{i \sqrt{b^2 c^2+b^2+1} (-1)^{\left\lfloor \frac{j-1}{2}\right\rfloor }+i b c (-1)^j}{\sqrt{b^2+1}},-\frac{b (-1)^i+i (-1)^{\left\lceil \frac{i-1}{2}\right\rceil }}{\sqrt{b^2+1}}} % \right. \\ &&\left. \right)\quad (ii) \end{eqnarray}$$

(*Definitions *)

Clear[F]; Clear[FF];
F[z_, a_, b_] :=
Log[a + z] Log[(b + z)/(-a + b)] + PolyLog[2, (a + z)/(a - b)];
FF[A_, B_, a_, b_] :=
Module[{result, ts, zs, zsp, zsm, eps = 10^(-30)},
(*This is Integrate[Log[z+a]/(z+b),{z,A,B}] where all a,b,A,
and B are complex. *)
result = F[B, a, b] - F[A, a, b];

ts = - (Im[(A + b) (Conjugate[b] - Conjugate[a])]/
Im[(B - A) (Conjugate[b] - Conjugate[a])]);
If[0 <= ts <= 1,
zsp = A + (ts + eps) (B - A);
zsm = A + (ts - eps) (B - A);
result += -F[zsp, a, b] + F[zsm, a, b];
];

result
];

rho[xi_] := Exp[-xi^2/2]/Sqrt[2 Pi]; a =.; b =.; c =.; eps = 10^(-12);
J[a_, b_, c_] :=
NIntegrate[rho[xi] Erf[a xi] OwenT[ b xi, c], {xi, 0, Infinity},
WorkingPrecision -> 20];

For[count = 1, count <= 100, count++,
{a, b, c} = RandomReal[{-10, 10}, 3, WorkingPrecision -> 50];
X1 = J[a, b, c];

X2 = 1/
Pi^2 (ArcTan[Sqrt[2] a]/2 ArcTan[ c] -
1/8  Sum[
FF[1, ( Sqrt[1 + 2 a^2 + b^2] - Sqrt[2] a)/Sqrt[
1 + b^2], ((-1)^j I b c + (-1)^Floor[(j - 1)/2] I Sqrt[
1 + b^2 + b^2 c^2])/Sqrt[
1 + b^2], -(((-1)^Ceiling[(i - 1)/2] I + (-1)^i b)/Sqrt[
1 + b^2])] (-1)^(j - 1 + Floor[(i - 1)/2]), {i, 1, 4}, {j,
1, 4}] );
If[Abs[X1/X2 - 1] > 10^(-6), Print[{a, b, c, X1, X2}]; Break[]];
If[Mod[count, 10] == 0, PrintTemporary[count]];
];
Print["All matches."];

All matches.


Now my question is what is the result for bigger values of $$\nu_1$$ and $$\nu_2$$ ?

Here we give an answer in the case $$\nu_1=2$$ and $$\nu_2=1$$. Firstly we have: $$\begin{eqnarray} \partial_{a_2} {\mathfrak J}^{(2,1)}(\vec{a},b,c) = \frac{a_1 \arcsin\left(\frac{c \sqrt{2 a_1^2+2 a_2^2+1}}{\sqrt{c^2+1} \sqrt{2 a_1^2+2 a_2^2+b^2+1}}\right)}{\pi ^2 \left(2 a_2^2+1\right) \sqrt{2 a_1^2+2 a_2^2+1}}-\frac{b \arctan\left(\frac{\sqrt{2} a_1 b c}{\sqrt{2 a_2^2+b^2+1} \sqrt{2 a_1^2+2 a_2^2+b^2 c^2+b^2+1}}\right)}{\sqrt{2} \pi ^2 \left(2 a_2^2+1\right) \sqrt{2 a_2^2+b^2+1}} \quad (i) \end{eqnarray}$$ The result comes from the definition of the error function, then from integrating by parts once with respect to $$\xi$$ and then using integral identities given in the two links in the formulation of the question above. Now since $${\mathfrak J}^{(2,1)}(a_1,0,b,c) = 0$$ all we need to do is to integrate the right hand side above with respect to $$a_2$$. I will show that all the integrals that emerge are feasible and are reduced to elementary function and to di-logarithms. We have: $$\begin{eqnarray} &&{\mathfrak J}^{(2,1)}(\vec{a},b,c) =\\ && \underbrace{\int\limits_0^{a_2} \frac{a_1 \arcsin\left(\frac{c \sqrt{2 a_1^2+2 \xi ^2+1}}{\sqrt{c^2+1} \sqrt{2 a_1^2+b^2+2 \xi ^2+1}}\right)}{\pi ^2 \left(2 \xi ^2+1\right) \sqrt{2 a_1^2+2 \xi ^2+1}}d\xi}_{{\mathfrak I}_1} -\underbrace{\int\limits_0^{a_2}\frac{b \arctan\left(\frac{\sqrt{2} a_1 b c}{\sqrt{b^2+2 \xi ^2+1} \sqrt{2 a_1^2+b^2 c^2+b^2+2 \xi ^2+1}}\right)}{\sqrt{2} \pi ^2 \left(2 \xi ^2+1\right) \sqrt{b^2+2 \xi ^2+1}}d\xi}_{{\mathfrak I}_2}=\\ && -\frac{\arctan\left(\frac{\sqrt{2} a_2 b}{\sqrt{2 a_2^2+b^2+1}}\right) \arctan\left(\frac{\sqrt{2} a_1 b c}{\sqrt{2 a_2^2+b^2+1} \sqrt{2 a_1^2+2 a_2^2+b^2 c^2+b^2+1}}\right)}{2 \pi ^2}+\\ &&\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\! sign(a_2) \frac{a_1 c^3 \imath \sqrt{1+c^2} |b|}{\pi^2 |c|} \int\limits_{\frac{\sqrt{1+2 a_1^2}c}{\sqrt{1+2 a_1^2+b^2(1+c^2)}}}^{\frac{\sqrt{1+2 a_1^2+2 a_2^2}c}{\sqrt{1+2 a_1^2+2 a_2^2+b^2(1+c^2)}}} \frac{\arctan(u)}{\sqrt{2 \left(2 a_1^2+1\right) c^2-2 u^2 \left(2 a_1^2+b^2 \left(c^2+1\right)+1\right)} \left(2 a_1^2 \left(u^2-c^2\right)+b^2 \left(c^2+1\right) u^2\right)} du+\\ &&\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!-\frac{\sqrt{2} ba_1 c \sqrt{1+b^2}}{2 \pi^2 |b|} \int\limits_{0}^{\frac{\sqrt{2} b a_2}{\sqrt{1+b^2+2 a_2^2}}} \frac{u \arctan(u) \left(b^2 \left(2 a_1^2-c^2 u^2+2\right)-2 a_1^2 u^2+b^4 \left(c^2+2\right)\right)}{\left(\left(2 a_1^2+1\right) b^2-2 a_1^2 u^2+b^4\right) \left(b^2 \left(c^2+1\right)-c^2 u^2+1\right) \sqrt{u^2 \left(-2 a_1^2-b^2 c^2\right)+b^2 \left(2 a_1^2+b^2 \left(c^2+1\right)+1\right)}} du % \end{eqnarray}$$ now, we obtained the second line in the following way. In the integral $${\mathfrak I}_1$$ we just substituted for the argument of the $$\arcsin$$ and then we used the identity $$\arcsin(x)= \arctan(x/\sqrt{1-x^2})$$ and again we substituted for the argument of the $$\arctan$$. In the integral $${\mathfrak I}_2$$ we integrated by parts using the identity : $$$$\int\frac{1}{(1+2 \xi^2) \sqrt{1+2 \xi^2+b^2}} d\xi = \frac{1}{\sqrt{2} b} \arctan\left( \frac{\sqrt{2} b \xi}{\sqrt{1+b^2+2 \xi^2}}\right)$$$$ and then we just substituted for the argument of the $$\arctan$$.

Now both integrals over $$u$$ are treated in the same way. Firstly we use a trigonometric substitution $$u= C \sin(\phi)$$ where $$C$$ is chosen in an appropriate way (meaning such that the square cancels the cosine term from the differential) and after that we use the well known substitution for $$\tan(\phi/2)$$. Finally we also use the identity $$\arctan(x) = 1/(2 \imath) \log((1+\imath x)/(1-\imath x))$$.

The final result is as follows: $$\begin{eqnarray} &&{\mathfrak J}^{(2,1)}(\vec{a},b,c)=\\ &&-\frac{1}{2 \pi^2} \arctan\left( \frac{\sqrt{2} b a_2}{\sqrt{1+b^2+2 a_2^2}}\right) \arctan\left( \frac{a_1 b c \sqrt{2}}{\sqrt{1+2 a_2^2+b^2} d_2}\right)+\\ &&-sign(a_2 c)\frac{\sqrt{1+c^2} |b|}{\sqrt{2} \pi^2 2 a_1 d_1} \int\limits_{sign(c)}^{\frac{1-\sqrt{1-x_1^2}}{x_1}} \frac{1+u^2}{(u^2-\frac{b\sqrt{2(1+c^2)}}{a_1 d_1} u-1)(u^2+\frac{b\sqrt{2(1+c^2)}}{a_1 d_1} u-1)} \cdot \log\left[ \frac{(u+c_1)(u+c_2)}{(u+c_3)(u+c_4)}\right]du+\\ &&\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\frac{i \sqrt{2} a_1 b \sqrt{b^2+1} c d_1 \left(2 a_1^2+b^2 c^2\right)}{\pi ^2 \left(2 a_1^2 d_1^2+b^2 \left(b^2+1\right) c^2\right) \left(2 a_1^2 \left(b^2+1\right)+b^2 c^2 d_1^2\right)}\cdot \int\limits_0^{\frac{1-\sqrt{1-x_2^2}}{x_2}} \frac{(2+2 a_1^2+b^2(2+c^2))u - 2(a_1^2+b^2 c^2) u^3 + (2+2 a_1^2 + b^2(2+c^2))u^5)}{\prod\limits_{j=1}^4 (u^2-r_j)} \cdot \log\left[ \frac{(u+c^{(1)}_1)(u+c^{(1)}_2)}{(u+c^{(1)}_3)(u+c^{(1)}_4)}\right]du \quad (iii) \end{eqnarray}$$ where: $$\begin{eqnarray} c&:=&\imath \left( C-\sqrt{1+C^2},C+\sqrt{1+C^2},-C+\sqrt{1+C^2},-C-\sqrt{1+C^2}\right)\\ c^{(1)}&:=&\imath \left( C^{(1)}-\sqrt{1+[C^{(1)}]^2},C^{(1)}+\sqrt{1+[C^{(1)}]^2},-C^{(1)}+\sqrt{1+[C^{(1)}]^2},-C^{(1)}-\sqrt{1+[C^{(1)}]^2}\right)\\ x_1&:=&sign(c) \frac{d_1 \sqrt{1+2 a_1^2+2 a_2^2}}{d_2 \sqrt{1+2 a_1^2}}\\ x_2&:=& sign(b) \frac{\sqrt{2 a_1^2+b^2 c^2}\sqrt{2} a_2}{d_1 \sqrt{1+b^2+2 a_2^2}}\\ r&:=& \left( \frac{(\imath b c\sqrt{1+b^2} - \sqrt{2} a_1 d_1)^2}{2 a_1^2 d_1^2 + b^2 c^2(1+b^2)},\frac{(-\imath b c\sqrt{1+b^2} - \sqrt{2} a_1 d_1)^2}{2 a_1^2 d_1^2 + b^2 c^2(1+b^2)},\frac{(\imath a_1 \sqrt{2(1+b^2)} - b c d_1)^2}{b^2 c^2 d_1^2 + 2 a_1^2(1+b^2)},\frac{(-\imath a_1 \sqrt{2(1+b^2)} - b c d_1)^2}{b^2 c^2 d_1^2 + 2 a_1^2(1+b^2)}\right) \end{eqnarray}$$ where $$\begin{eqnarray} C&:=& \frac{\sqrt{(1+2 a_1^2) c^2}}{d_1}\\ C^{(1)}&:=& \frac{|b| d_1}{\sqrt{2 a_1^2+b^2 c^2}} \end{eqnarray}$$ and $$\begin{eqnarray} d_1&:=& \sqrt{1+2 a_1^2+b^2(1+c^2)}\\ d_2&:=& \sqrt{1+2 a_1^2+2 a_2^2+b^2(1+c^2)} \end{eqnarray}$$

Now it is clear that the integrals over $$u$$ in $$(iii)$$ can always be reduced to di-logarithms by decomposing the rational function in the integrand into partial fractions and then using the product properties of the logarithms. We could have formally carried out this task but this would have made the expressions even more unwieldy and would not bring much insight at all. Yet having said this we can say that the computation is completed. As usual I enclose a piece of code that verifies numerically the calculations above.

ll = {};
For[count = 1, count <= 100, count++,
{a1, a2, b, c} = RandomReal[{-1, 1}, 4, WorkingPrecision -> 50];
I1 = NIntegrate[
rho[xi]  Erf[a1 xi] Erf[a2 xi] OwenT[b xi, c], {xi, 0, Infinity},
WorkingPrecision -> 20];
I2 = NIntegrate[((
2 a1 ArcSin[(Sqrt[1 + 2 a1^2 + 2 xi^2] c)/(
Sqrt[1 + 2 a1^2 + 2 xi^2 + b^2] Sqrt[1 + c^2])])/Sqrt[
1 + 2 a1^2 + 2 xi^2] - (
Sqrt[2] b ArcTan[(a1 b c Sqrt[2])/(
Sqrt[(1 + 2 xi^2 + b^2)] Sqrt[
1 + 2 a1^2 + 2 xi^2 + b^2 + b^2 c^2])])/Sqrt[
1 + 2 xi^2 + b^2])/(2 (1 + 2 xi^2) \[Pi]^2), {xi, 0, a2}];
I3 = Sign[a2] NIntegrate[(
a1 c^3 Sqrt[-1 - c^2] Abs[b] ArcTan[v])/(\[Pi]^2 Sqrt[
2 (1 + 2 a1^2) c^4 -
2 c^2 (1 + 2 a1^2 + b^2 (1 + c^2)) v^2] (b^2 (1 + c^2) v^2 +
2 a1^2 (-c^2 + v^2))), {v, A1/Sqrt[1 - A1^2] , A2/Sqrt[
1 - A2^2] }] -
Sqrt[2] b/(2  \[Pi]^2) NIntegrate[
ArcTan[(a1 b c Sqrt[2])/(
Sqrt[(1 + 2 xi^2 + b^2)] Sqrt[
1 + 2 a1^2 + 2 xi^2 + b^2 + b^2 c^2])]/( (1 + 2 xi^2) Sqrt[
1 + 2 xi^2 + b^2]), {xi, 0, a2}];
I4 = -(1/(2  \[Pi]^2)) ArcTan[(Sqrt[2] b a2)/Sqrt[
1 + b^2 + 2 a2^2]] ArcTan[(a1 b c Sqrt[2])/(
Sqrt[(1 + 2 a2^2 + b^2)] Sqrt[
1 + 2 a1^2 + 2 a2^2 + b^2 + b^2 c^2])] +
Sign[a2] (a1 c^3 Sqrt[-1 - c^2] Abs[b])/(\[Pi]^2 Abs[c])
NIntegrate[
ArcTan[v]/((b^2 (1 + c^2) v^2 + 2 a1^2 (-c^2 + v^2)) Sqrt[
2 (1 + 2 a1^2) c^2 -
2 (1 + 2 a1^2 + b^2 (1 + c^2)) v^2] ), {v, (
Sqrt[1 + 2 a1^2] c)/ Sqrt[1 + 2 a1^2 + b^2 (1 + c^2)], (
Sqrt[1 + 2 a1^2 + 2 a2^2] c)/ Sqrt[
1 + 2 a1^2 + 2 a2^2 + b^2 + b^2 c^2]}] + -Sqrt[2]
b/(2  \[Pi]^2) (a1 Sqrt[1 + b^2] c)/
Abs[b] NIntegrate[(
u  (b^4 (2 + c^2) - 2 a1^2 u^2 +
b^2 (2 + 2 a1^2 - c^2 u^2)) ArcTan[
u])/( ((1 + 2 a1^2) b^2 + b^4 - 2 a1^2 u^2) (1 +
b^2 (1 + c^2) - c^2 u^2) Sqrt[
b^2 (1 + 2 a1^2 + b^2 (1 + c^2)) + (-2 a1^2 -
b^2 c^2) u^2]), {u, 0, (Sqrt[2] b a2)/Sqrt[
1 + b^2 + 2 a2^2]}];
(*Now a trigonometric substitution u--> CC Sin[
phi] followed by the u = Tan[phi/2] substitution. *)
{d1, d2} = {Sqrt[1 + 2 a1^2 + b^2 (1 + c^2)], Sqrt[
1 + 2 a1^2 + 2 a2^2 + b^2 + b^2 c^2]};
{CC, CC1} = { Sqrt[(1 + 2 a1^2) c^2]/ d1, (Abs[b] d1)/
Sqrt[(2 a1^2 + b^2 c^2)]};
cc = {I (CC - Sqrt[1 + CC^2]), I (CC + Sqrt[1 + CC^2]),
I (-CC + Sqrt[1 + CC^2]), I (-CC - Sqrt[1 + CC^2])};
cc1 = {I (CC1 - Sqrt[1 + CC1^2]), I (CC1 + Sqrt[1 + CC1^2]),
I (-CC1 + Sqrt[1 + CC1^2]), I (-CC1 - Sqrt[1 + CC1^2])};
{x1, x2} = {Sign[c] (d1 Sqrt[1 + 2 a1^2 + 2 a2^2] )/(
d2 Sqrt[ (1 + 2 a1^2)]),
Sign[b] (Sqrt[(2 a1^2 + b^2 c^2)] Sqrt[2] a2)/(
d1 Sqrt[1 + b^2 + 2 a2^2])};
rr = {(I b c Sqrt[1 + b^2] - Sqrt[2] a1 d1)^2/(
2 a1^2 d1^2 +
b^2 c^2 (1 + b^2)), (-I b c Sqrt[1 + b^2] - Sqrt[2] a1 d1)^2/(
2 a1^2 d1^2 +
b^2 c^2 (1 + b^2)), (I a1 Sqrt[2 (1 + b^2)] - b c d1)^2/(
b^2 c^2 d1^2 +
2 (a1^2) (1 + b^2) ), (-I a1 Sqrt[2 (1 + b^2)] - b c d1)^2/(
b^2 c^2 d1^2 + 2 (a1^2) (1 + b^2) )};

I5 = -(1/(2  \[Pi]^2)) ArcTan[(Sqrt[2] b a2)/Sqrt[
1 + b^2 + 2 a2^2]] ArcTan[(a1 b c Sqrt[2])/(
Sqrt[(1 + 2 a2^2 + b^2)] d2)] + -Sign[a2 c] (
Sqrt[1 + c^2] Abs[b])/(Sqrt[2]  \[Pi]^2 2 a1 d1)
NIntegrate[(1 +
u^2)/((u^2 - (b Sqrt[2 (1 + c^2)])/(a1 d1) u -
1) (u^2 + (b Sqrt[2 (1 + c^2)])/(a1 d1) u - 1))
Log[((cc[[1]] + u) (cc[[2]] + u))/((cc[[3]] + u) (cc[[4]] +
u))], {u, Sign[c], (1 - Sqrt[1 - x1^2])/x1},
WorkingPrecision -> 20] + (
I Sqrt[2] a1 c b d1 Sqrt[
1 + b^2] (2 a1^2 + b^2 c^2))/(\[Pi]^2 (2 a1^2 d1^2 +
b^2 c^2 (1 + b^2)) (b^2 c^2 d1^2 + 2 a1^2 (1 + b^2)))
NIntegrate[((2 + 2 a1^2 + b^2 (2 + c^2)) u -
2 (2 a1^2 + b^2 c^2) u^3 + (2 + 2 a1^2 +
b^2 (2 + c^2)) u^5)/((u^2 - rr[[1]]) (u^2 - rr[[2]]) (u^2 -
rr[[3]]) (u^2 - rr[[4]]))
Log[((cc1[[1]] + u) (cc1[[2]] + u))/((cc1[[3]] +
u) (cc1[[4]] + u))], {u, 0, (1 - Sqrt[1 - x2^2])/x2},
WorkingPrecision -> 20];
If[Abs[I2/I1 - 1] > 10^(-3), Print[{count, {a1, a2, b, c, I1, I2}}];
Break[]];
If[Mod[count, 10] == 0, PrintTemporary[count]];
ll = Join[ll, {{I1, I2, I3, I4, I5}}];
];
Abs[ll[[All, 1]]/ll[[All, -1]] - 1]