An integral involving a Gaussian and an Owen's T function. Let $a,b,c \ge 0$.
In the course of solving gaussian integral of power of cdf : $\int_{-\infty}^{+\infty} \Phi(x)^n \cdot \phi(a+bx) \cdot dx$ we came across a following integral:
\begin{eqnarray}
T(a,b,c):= \int\limits_a^\infty \phi(\xi) T(\xi b,c) d\xi
\end{eqnarray}
where $T(.,.)$ is the Owen's T function https://en.wikipedia.org/wiki/Owen%27s_T_function. As the Wikipedia page tells us this function is related to likelihoods of a bivariate Gaussian distribution and furthermore the integral itself is related to likelihoods of a triple-variate Gaussian distribution.
Now, by going to spherical coordinates we managed to find our integral in case $a=0$. We have:
\begin{equation}
T(0,b,c) = \frac{1}{4\pi} \arcsin\left( 
\frac{c}{\sqrt{(1+b^2)(1+c^2)}}\right)
\end{equation}
My question is therefore what is the value of the integral for generic $a>0$ ?

It turns out that the answer is actually quite simple and can be obtained using elementary methods. We have:
\begin{eqnarray}
T(a,b,c)= \int\limits_{{\mathbb R}^3} 1_{\xi_0 > a} 1_{\xi_1 > \xi_0 b} 1_{c \xi_1 > \xi_2 >0} \prod\limits_{p=0}^2 \phi(\xi_p) d \xi_p
\end{eqnarray}
Now we go to spherical coordinates as follows:
\begin{eqnarray}
\xi_0&=&r \sin(\theta) \cos(\phi)\\
\xi_1&=&r \sin(\theta) \sin(\phi)\\
\xi_2&=& r \cos(\theta)
\end{eqnarray}
where $\theta \in[0,\pi/2)$ and $\phi\in[-\pi/2,\pi/2)$ which follows from the fact that all three coordinates are positive.
Now we carefuly analyze the inequality conditions:
\begin{eqnarray}
1_{\xi_0 > a} &:& r \sin(\theta) \cos(\phi) > a \quad \Longrightarrow \quad \theta > \frac{a}{r \cos(\phi)} \quad \Longrightarrow \quad \theta > \arcsin(\frac{a}{r \cos(\phi)})\\
1_{\xi_1>\xi_0 b} &:& r \sin(\theta) \sin(\phi) > b r \sin(\theta) \cos(\phi) \quad \Longrightarrow \quad \tan(\phi) > b\\
1_{c \xi_1> \xi_2 >0} &:& c r \sin(\theta) \sin(\phi) > r \cos(\theta) \quad \Longrightarrow \quad c \sin(\phi) > \cot(\theta) \quad \Longrightarrow \quad \theta > arccot(c \sin(\phi))
\end{eqnarray}
Now, in the last equation on the right in the top line above we must have $r > a/\cos(\phi) = a \sqrt{1+\tan(\phi)^2} > a \sqrt{1+b^2}$.
Therefore we write down the integral in spherical coordinates as follows:
\begin{eqnarray}
(2\pi)^{3/2}T(a,b,c)&=& \int\limits_{a\sqrt{1+b^2}}^\infty e^{-1/2 r^2} r^2 \int\limits_{\arctan(b)}^{\arccos(a/r)} \int\limits_{\arcsin(\frac{a}{(r \cos(\phi))})\vee arccot(c \sin(\phi))}^{\pi/2} \sin(\theta) d\theta \quad d\phi \quad dr\\
&=& \int\limits_{a\sqrt{1+b^2}}^\infty e^{-1/2 r^2} r^2 \int\limits_{\arctan(b)}^{\arccos(a/r)} 
%
\left[\cos(\arcsin(\frac{a}{r \cos(\phi)})) \wedge 
\cos(arccot(c \sin(\phi)))
\right]
%
\quad d\phi \quad dr\\
&=&
\int\limits_{a\sqrt{1+b^2}}^\infty e^{-1/2 r^2} r^2 \int\limits_{\arctan(b)}^{\arccos(a/r)} 
%
\left[
\frac{c \sin(\phi)}{\sqrt{1+c^2 \sin(\phi)^2}} \wedge
\frac{\sqrt{r^2 \cos(\phi)^2-a^2}}{r \cos(\phi)}
\right]
%
\quad d\phi \quad dr  \\
&=&
\int\limits_{a\sqrt{1+b^2}}^{a\sqrt{1+b^2(1+c^2)}} e^{-1/2 r^2} r^2 
\int\limits_{\arctan(b)}^{\arccos(a/r)}
\frac{\sqrt{r^2 \cos(\phi)^2-a^2}}{r \cos(\phi)} d\phi dr +\\
&&
\int\limits_{a\sqrt{1+b^2(1+c^2)}}^\infty e^{-1/2 r^2} r^2 
\left(
\int\limits_{\arctan(b)}^{\arccos(\frac{a \sqrt{1+c^2}}{\sqrt{a^2 c^2 + r^2}})}
\frac{c \sin(\phi)}{\sqrt{1+c^2 \sin(\phi)^2}}
d\phi + 
%
\int\limits_{\arccos(\frac{a \sqrt{1+c^2}}{\sqrt{a^2 c^2 + r^2}})}^{\arccos(a/r)}
\frac{\sqrt{r^2 \cos(\phi)^2-a^2}}{r \cos(\phi)} d\phi
\right) =\\
&=&
\int\limits_{a\sqrt{1+b^2}}^{a\sqrt{1+b^2(1+c^2)}} e^{-1/2 r^2} r^2
\left(
\frac{a}{r} \arctan[\frac{a b}{\sqrt{r^2-a^2(1+b^2)}}]-
\arctan[\frac{r b}{\sqrt{r^2-a^2(1+b^2)}}]
\right)dr +\\
&& \frac{\pi ^{3/2}}{2 \sqrt{2}} \left(\sqrt{\frac{2}{\pi }} a \left(\sqrt{b^2+1}-1\right) e^{-\frac{1}{2} a^2 \left(b^2+1\right)}-\frac{\sqrt{2} a e^{-\frac{1}{2} a^2 \Delta^2} \left(\Delta \left(\pi -2 \arctan\left(\frac{c}{\Delta}\right)\right)-2 \arctan\left(\frac{1}{c}\right)\right)}{\pi ^{3/2}}-\text{erf}\left(\frac{a \sqrt{b^2+1}}{\sqrt{2}}\right)+\frac{\text{erf}\left(\frac{a \Delta}{\sqrt{2}}\right)
   \left(\pi -2 \arctan\left(\frac{c}{\Delta}\right)\right)}{\pi }+\frac{2 \arctan\left(\frac{c}{\Delta}\right)}{\pi }\right)
\end{eqnarray}
where $\Delta:=\sqrt{1+b^2(1+c^2)}$.
As a sanity check we look at the case $a=0$. In here only the very last term  survives and we have:
\begin{eqnarray}
(2\pi)^{3/2} T(0,b,c)&=& \sqrt{\frac{\pi}{2}} \arctan(\frac{c}{\Delta})\\
\Longrightarrow\\
T(0,b,c)&=& \frac{1}{4\pi} \arctan(\frac{c}{\Delta})= \frac{1}{4\pi} \arcsin(\frac{c}{\sqrt{(1+b^2)(1+c^2)}})
\end{eqnarray}
as it should be.
