As Wikipedia teaches us https://en.wikipedia.org/wiki/Owen%27s_T_function the Owen's T function $T(h,a)$ defines a probability of a bivariate event $X>h$ and $0<Y<a X$ where $X,Y$ are standard, independent Gaussian random variables.
Now in the context of question Multivariate gaussian integral over positive reals a necessity appeared to deal with a slightly more general quantity. \begin{equation} T(h,a,b):= {\bf P}\left(X>h \quad \wedge \quad a X+b > Y > 0 \left. \right| X = N(0,1) , Y=N(0,1) \right) \end{equation} We have shown that : \begin{eqnarray} &&T(h,a,b)= \int\limits_h^\infty \frac{\exp(-1/2 \xi^2)}{\sqrt{2\pi}} \frac{1}{2} Erf(\frac{a \xi+b}{\sqrt{2}}) d\xi \quad (i1)\\ &&= \int\limits_0^a \frac{e^{-\frac{b^2}{2}-b h \xi -\frac{1}{2} h^2 \left(\xi ^2+1\right)}}{2 \pi \left(\xi ^2+1\right)} d\xi - \frac{b}{2\sqrt{2}\sqrt{\pi}} \int\limits_0^a \frac{\xi e^{-\frac{b^2}{2 \xi ^2+2}} \text{erfc}\left(\frac{\xi (b+h \xi )+h}{\sqrt{2} \sqrt{\xi ^2+1}}\right)}{\left(\xi ^2+1\right)^{3/2}} d\xi + \frac{1}{4} \text{erf}\left(\frac{b}{\sqrt{2}}\right) \text{erfc}\left(\frac{h}{\sqrt{2}}\right) \quad (i2) \end{eqnarray}
{a, b, h} = RandomReal[{0, 1}, 3, WorkingPrecision -> 50]; b = 0;
NIntegrate[
Exp[-x^2/2]/Sqrt[2 Pi] 1/2 Erf[(a x + b)/Sqrt[2]], {x, h, Infinity},
WorkingPrecision -> 20]
NIntegrate[(E^(-(b^2/2) - xi b h - 1/2 (1 + xi^2) h^2)) /(
2 (1 + xi^2) \[Pi]) -
b /(2 Sqrt[2] Sqrt[ \[Pi]]) (
xi Erfc[(h + xi (b + xi h))/(Sqrt[2] Sqrt[1 + xi^2])])/ ((1 +
xi^2)^(3/2)) E^(-(b^2/(2 + 2 xi^2))), {xi, 0, a},
WorkingPrecision -> 20] + Erfc[h/Sqrt[2]] Erf[b/Sqrt[2]] 1/4
Update: Let $A_j \in {\mathbb R}$ for $j=1,\cdots,3$ and let $x\in {\mathbb R}$. Then we have: \begin{eqnarray} T(A_1 x, A_2, A_3 x) = \frac{1}{2\pi} \left(\arctan(A_2)-\arctan(A_2+\frac{A_3}{A_1})-\arctan(\frac{A_1+A_2 A_3+A_2^2 A_1}{A_3})\right) + \frac{1}{4} erf[\frac{A_3 x}{\sqrt{2} \sqrt{1+A_2^2}}] + T(A_1 x, \frac{A_2 A_1+A_3}{A_1})+ T(\frac{A_3 x}{\sqrt{1+A_2^2}},\frac{A_1+A_2 A_3 + A_2^2 A_1}{A_3}) \quad (ii) \end{eqnarray}
This identity comes from differentiating both sides with respect to $x$ then using the definition of the generalized Owen's T function to evaluate the derivative on the right hand side and having done this integrating both sides with respect to $x$ again.
Let us present the proof of that in detail. Firstly we define $f(x) := T[A_1 x, A_2, A_3 x]$. Now we compute the derivative using the chain rule. We have: \begin{eqnarray} \frac{d }{d x} f(x) &=& \partial_1 T[A_1 x, A_2 , A_3 x] \cdot A_1 + \partial_3 T[A_1 x, A_2, A_3 x] \cdot A_3 \\ &=& - \left. \rho(h) \frac{1}{2} erf[\frac{a h+b}{\sqrt{2}}] \right|_{\begin{array}{r} h=A_1 x \\ a=A_2 \\ b=A_3 x \end{array}}\cdot A_1 + \left.\frac{1}{\sqrt{1+a^2}} \frac{1}{2} erf[\frac{h+a^2 h+a b}{\sqrt{2} \sqrt{1+a^2}}] \rho(\frac{b}{1+a^2})\right|_{\begin{array}{r} h=A_1 x \\ a=A_2 \\ b=A_3 x \end{array}} \cdot A_3 \\ &=& -\rho(A_1 x) \frac{1}{2} erf[\frac{A_1 A_2 + A_3}{\sqrt{2}} x] \cdot A_1 + \frac{1}{\sqrt{1+A_2^2}} \rho(\frac{A_3 x}{\sqrt{1+A_2^2}}) \frac{1}{2} erfc[\frac{A_1+A_2 A_3 +A_1 A_2^2}{\sqrt{2} \sqrt{1+A_2^2}} x] \cdot A_3 \end{eqnarray}
Now we integrate. We have: \begin{eqnarray} f(x)- f(0) &=& - \int\limits_0^x \rho(A_1 \xi) \frac{1}{2} erf[\frac{A_1 A_2 + A_3}{\sqrt{2}} \xi] d\xi \cdot A_1 + \\ &&\frac{1}{\sqrt{1+A_2^2}}\int\limits_0^x \rho(\frac{A_3 \xi}{\sqrt{1+A_2^2}}) \frac{1}{2} erfc[\frac{A_1+A_2 A_3 +A_1 A_2^2}{\sqrt{2} \sqrt{1+A_2^2}} \xi] d\xi \cdot A_3 \\ f(x) - \frac{1}{2\pi} \arctan(A_2) &=& - \frac{1}{2\pi} \arctan\left( \frac{A_1 A_2+A_3}{A_1}\right) + T(A_1 x, \frac{A_1 A_2 + A_3}{A_1}) + \\ && \frac{1}{4} erf\left( \frac{A_3}{\sqrt{2} \sqrt{1+A_2^2}} x\right) +\\ &&-\frac{1}{2\pi} \arctan\left( \frac{A_1+A_2 A_3 + A_1 A_2^2}{A_3}\right) + T\left( \frac{A_3}{\sqrt{1+A_2^2}} x, \frac{A_1+A_2 A_3 + A_1 A_2^2}{A_3}\right) \end{eqnarray} where in the second line we used the results from An integral involving error functions and a Gaussian and the definition of the Owen's T function. This completes the proof.
(*A certain derivative. Used in Q869502.nb*)
T[h_, a_, b_] :=
NIntegrate[(E^(-(b^2/2) - xi b h - 1/2 (1 + xi^2) h^2)) /(
2 (1 + xi^2) \[Pi]) -
b /(2 Sqrt[2] Sqrt[ \[Pi]]) (
xi Erfc[(h + xi (b + xi h))/(Sqrt[2] Sqrt[1 + xi^2])])/ ((1 +
xi^2)^(3/2)) E^(-(b^2/(2 + 2 xi^2))), {xi, 0, a},
WorkingPrecision -> 20] + Erfc[h/Sqrt[2]] Erf[b/Sqrt[2]] 1/4;
{A1, A2, A3} = RandomReal[{-1, 1}, 3, WorkingPrecision -> 50];
u = Range[0, 1, 1/100];
mT = Interpolation[Transpose[{u, T[A1 u, A2, A3 u]}]];
u =.; u = RandomReal[{0, 1}, WorkingPrecision -> 50];
mT'[u]
-rho[A1 u] 1/2 Erf[(A1 A2 + A3)/Sqrt[2] u] A1 +
1/Sqrt[1 + A2^2]
rho[(A3 u)/Sqrt[1 + A2^2]] 1/
2 Erfc[(A1 + A2 A3 + A1 A2^2)/(Sqrt[2] Sqrt[1 + A2^2]) u] A3
T[A1 u, A2, A3 u]
1/(2 Pi) (ArcTan[A2] - ArcTan[(A2 A1 + A3)/A1] -
ArcTan[(A1 + A2 A3 + A2^2 A1)/A3]) +
1/4 Erf[(A3 u)/(Sqrt[2] Sqrt[1 + A2^2])] +
OwenT[A1 u, (A2 A1 + A3)/A1] +
OwenT[A3/Sqrt[1 + A2^2] u, (A1 + A2 A3 + A2^2 A1)/A3]
1/(2 Pi) (-ArcTan[A3/((A1 + A2 A3 + A2^2 A1))] -
ArcTan[(A1 + A2 A3 + A2^2 A1)/A3]) +
1/4 Erf[(A3 u)/(Sqrt[2] Sqrt[1 + A2^2])] +
OwenT[A1 u, (A2 A1 + A3)/A1] +
OwenT[A3/Sqrt[1 + A2^2] u, (A1 + A2 A3 + A2^2 A1)/A3]
-1/(2 Pi) Pi/2 (Sign[A3/((A1 + A2 A3 + A2^2 A1))]) +
1/4 Erf[(A3 u)/(Sqrt[2] Sqrt[1 + A2^2])] +
OwenT[A1 u, (A2 A1 + A3)/A1] +
OwenT[A3/Sqrt[1 + A2^2] u, (A1 + A2 A3 + A2^2 A1)/A3]
-(1/4) Sign[A3/((A1 + A2 A3 + A2^2 A1))] +
1/4 Erf[(A3 u)/(Sqrt[2] Sqrt[1 + A2^2])] +
OwenT[A1 u, (A2 A1 + A3)/A1] +
OwenT[A3/Sqrt[1 + A2^2] u, (A1 + A2 A3 + A2^2 A1)/A3]
Now by both taking $x=1$ and replacing $A_1$,$A_2$ and $A_3$ by $h$, $a$ and $b$ in $(ii)$ we express the generalized Owen's T function through Owen's T function itself. We have: \begin{eqnarray} T(h,a,b) = \frac{1}{2\pi} \left(\arctan(a)-\arctan(a+\frac{b}{h})-\arctan(\frac{h+a b+a^2 h}{b})\right) + \frac{1}{4} erf[\frac{b}{\sqrt{2(1+a^2)}}] + T\left( h,\frac{a h+b}{h}\right) + T\left( \frac{b}{\sqrt{1+a^2}},\frac{h+a b+a^2 h}{b}\right) \end{eqnarray}
As a sanity check we look at the limit $b$ going to zero. We have: \begin{eqnarray} \lim_{b \rightarrow 0_+} T(h,a,b) &=& \frac{1}{2\pi} \left(\arctan(a)-\arctan(a)-\frac{\pi}{2} sign(h))\right) + 0 + T(h,a) + \frac{1}{4} sign(h) \\ &=& T(h,a) \end{eqnarray} as it should be.
As another sanity check we look at the case $a=\imath$. Going back to the calculations of the derivative above we have: \begin{eqnarray} \frac{d}{d x} f(x)= -\phi(A_1 x) \frac{1}{2} erf(\frac{A_1 A_2+A_3}{\sqrt{2}} x) A_1 + \frac{1}{2\pi \imath x} \exp(-\frac{1}{2} x^2 (2 A_1 \imath A_3+A_3^2)) \end{eqnarray} where we used the asymptotic expansion for the complementary error function given in https://en.wikipedia.org/wiki/Error_function#Complementary_error_function . Now we take a number $M$ such that $1< M$ and we integrate the above from unity to $M$ and we get: \begin{eqnarray} f(1)-f(M)= \left.\left( T(A_1 \cdot \xi,A_2+\frac{A_3}{A_1}) + \frac{1}{4\pi \imath} Ei(-\frac{1}{2}(1+2\imath \frac{A_1}{A_3})(\xi A_3)^2\right)\right|_{\xi=M}^{\xi=1} \end{eqnarray} where $Ei()$ is the exponential integral. Now it turns out that as $M\rightarrow \infty$ both $f(M)$ and $T(\dots M,\dots)$ tend to zero and \begin{equation} \lim\limits_{M\rightarrow \infty} \frac{1}{4 \pi \imath} Ei((a+\imath b) M) = sign(b) \cdot \frac{1}{4} \cdot 1_{a<0} + \infty \cdot 1_{a>0} \end{equation} Defining $b:=b_1+\imath b_2$ and taking $h>0$ this gives the final result: \begin{eqnarray} &&T(h,\imath , b) = \\ &&\left\{ \begin{array}{rr} T(h,\imath + \frac{b}{h}) + \frac{1}{4\pi \imath} Ei(\frac{1}{2}(-b_1^2+b_2^2+2 b_2 h-2\imath b_1(b_2+h))) + sign(b_1(b_2+h)) \cdot \frac{1}{4} & \mbox{if $b_2<0$ and $-b_1^2 + b_2^2+2 b_2 h <0$} \\ \infty & \mbox{otherwise} \end{array} \right. \end{eqnarray}
My question is the following. Has this quantity been ever analyzed in the literature before?
This is not an answer to the question above but instead it is a generalization of the results above. Define $\vec{a}:=(a_j)_{j=1}^d \in {\mathbb R}_+^d$ and let us define as $T^{(d)}(h,\vec{a})$ the probability of a following multivariate event $X>h$ and $0< Y_j < a_j X$ for $j=1,\cdots,d$ where $X$ and $\left( Y_j \right)_{j=1}^d$ are standard, independent Gaussian random variables.
Now take another vector $\vec{b}:=(b_j)_{j=1}^d \in {\mathbb R}_+^d$ and define a slightly more general quantity: \begin{eqnarray} T^{(d)}(h,\vec{a},\vec{b}) &:=& P\left( X > h \quad \wedge \quad \begin{array}{rrr} a_1 X + b_1 >& Y_1& >0 \\ \vdots &\vdots& \vdots \\ a_d X + b_d >& Y_d& >0 \end{array}\left. \right| \begin{array}{rrr} X&=&N(0,1)\\Y_1 &=&N(0,1)\\ &\vdots& \\ Y_d&=&N(0,1) \end{array} \right)\\ &=& \int\limits_h^\infty \rho(\xi) \left[\prod\limits_{i=1}^d \frac{1}{2} erf(\frac{a_i \xi+b_i}{\sqrt{2}})\right] d\xi \end{eqnarray} In what follows we will prove that if $d\le 2$ then the quantity $T^{(d)}(h,\vec{a},\vec{b})$ reduces to elementary functions and to $T^{(d)}(h,\vec{a})$ only.
As in the question above we consider a following quantity $T^{(d)}(h \cdot x, \vec{a}, \vec{b} \cdot x)$ which we differentiate with respect to $x$. We have: \begin{eqnarray} &&\frac{d}{d x} T^{(d)}(h \cdot x, \vec{a}, \vec{b} \cdot x)=\\ &&-h \cdot \rho(h \cdot x) \prod\limits_{i=1}^d \frac{1}{2} erf(\frac{a_i h x+b_i x}{\sqrt{2}})+\\ &&\sum\limits_{i=1}^d \frac{b_i}{\sqrt{2 \pi}} \int\limits_{h \cdot x}^\infty e^{-\frac{1}{2} (a_i \xi + b_i x)^2}\rho(\xi) \left[\prod\limits_{j=1,j\neq i}^d \frac{1}{2} erf(\frac{a_j h x+b_j x}{\sqrt{2}})\right] d\xi \end{eqnarray} What we do now is to simplify the second term on the right hand side, i.e. absorb the exponential into the Gaussian density and extract a constant prefactor. After that we integrate the identity above over $x$ from zero to unity. The result reads: \begin{eqnarray} &&T^{(d)}(h,\vec{a},\vec{b}) - T^{(d)}(0,\vec{a},\vec{0})=\\ &&-T^{(d)}(0,\vec{a}+\frac{1}{h} \vec{b},\vec{0}) + T^{(d)}(h,\vec{a}+\frac{1}{h} \vec{b},\vec{0}) + \\ &&\sum\limits_{i=1}^d \int\limits_0^{\frac{b_i}{\sqrt{1+a_i^2}}} \rho(x) \cdot T^{(d-1)}\left([a_i + \frac{h}{b_i}(1+a_i^2)]x,\frac{(a_j)_{j=1,j\neq i}^d}{\sqrt{1+a_i^2}},\frac{(b_j(1+a_i^2)-b_i a_i a_j)_{j=1,j\neq i}^d}{b_i \sqrt{1+a_i^2}} \right) dx \end{eqnarray} this clearly gave us a recurrence relation for the quantity in question subject to $T^{(0)}(h,\vec{a},\vec{b})=T^{(0)}(h)= 1/2 erfc(h/\sqrt{2})$.
Now we state the result for $d=2$. Firstly we define auxiliary quantities: \begin{eqnarray} \delta&:=& h^2+(a_1 h+b_1)^2+(a_2 h+b_2)^2\\ \delta_1&:=&h(1+a_1^2+a_2^2) + a_1 b_1 + a_2 b_2\\ \delta_2&:=&1+a_1^2+a_2^2\\ \hline\\ (m_1,m_2)&:=& (b_1(1+a_2^2)-a_1 a_2 b_2,b_2(1+a_1^2)-a_1 a_2 b_1)\\ (n_1,n_2)&:=&(h+h a_1^2+a_1 b_1,h+h a_2^2+a_2 b_2)\\ (o_1,o_2)&:=&(h a_1+b_1,h a_2 +b_2)\\ (p_1,p_2)&:=&(1+a_1^2,1+a_2^2)\\ (k_1,k_2)&:=&(\frac{\sqrt{p_1} \delta_1}{m_2},\frac{\sqrt{p_2} \delta_1}{m_1})\\ (l_1,l_2)&:=&(\frac{m_1}{\sqrt{p_2 \delta_2}},\frac{m_2}{\sqrt{p_1 \delta_2}}) \end{eqnarray} Then the result reads: \begin{eqnarray} &&4 \pi T^{(2)}(h,\vec{a},\vec{b})=\\ &&\arctan(\frac{a_1 a_2}{\sqrt{\delta_2}}) - \arctan(\frac{o_1 o_2}{h \sqrt{\delta}})+\\ &&\arctan(\frac{m_2}{\sqrt{\delta_2} b_1}) + \arctan(\frac{m_1}{\sqrt{\delta_2} b_2})+\\ &&\arctan(\frac{b_1 o_2}{n_1 \sqrt{\delta}}) + \arctan(\frac{b_2 o_1}{n_2 \sqrt{\delta}})+\\ &&\arctan(\frac{b_1 \delta_1}{\sqrt{m_2^2 \delta}}) + \arctan(\frac{b_2 \delta_1}{\sqrt{m_1^2 \delta}})+\\ &&\left( \arctan(\frac{a_2}{\sqrt{p_1}}) - \arctan(\frac{\sqrt{p_1} o_2}{n_1})-\arctan(k_1)\right) \cdot erf(\frac{b_1}{\sqrt{2 p_1}})+\\ &&\left( \arctan(\frac{a_1}{\sqrt{p_2}}) - \arctan(\frac{\sqrt{p_2} o_1}{n_2})-\arctan(k_2)\right) \cdot erf(\frac{b_2}{\sqrt{2 p_2}})+\\ &&2\pi \left( T(\frac{n_1}{\sqrt{p_1}},\frac{\sqrt{p_1} o_2}{n_1}) + T(l_2,k_1)\right)\cdot erf(\frac{b_1}{\sqrt{2 p_1}})+\\ &&2\pi \left( T(\frac{n_2}{\sqrt{p_2}},\frac{\sqrt{p_2} o_1}{n_2}) + T(l_1,k_2)\right)\cdot erf(\frac{b_2}{\sqrt{2 p_2}})+\\ &&-2\pi\left( T(\frac{b_1}{\sqrt{p_1}},\frac{m_2}{\sqrt{\delta_2} b_1}) + T(\frac{b_2}{\sqrt{p_2}},\frac{m_1}{\sqrt{\delta_2} b_2})\right)+\\ &&4\pi \left( T^{(2)}(h,(a_j+\frac{b_j}{h})_{j=1}^2) -T^{(2)}(\frac{n_1}{\sqrt{p_1}},(\frac{b_1}{n_1},\frac{\sqrt{p_1} o_2}{n_1})) -T^{(2)}(\frac{n_2}{\sqrt{p_2}},(\frac{b_2}{n_2},\frac{\sqrt{p_2} o_1}{n_2})) -T^{(2)}(l_2,(\frac{\sqrt{\delta_2}b_1}{m_2},k_1)) -T^{(2)}(l_1,(\frac{\sqrt{\delta_2}b_2}{m_1},k_2)) \right) \end{eqnarray} As usual I include a piece of code that verifies this expression:
Update: It might be interesting to know if it is possible to express the quantities $T^{(2)}(h,(a_1,a_2))$ in some alternative way. As a matter of fact by starting from the integral definition of this quantity , then differentiating with respect to $a_1$ and then integrating by parts and finally integrating with respect to $a_1$ from zero to $a_1$ we stumbled on a following formula: \begin{eqnarray} T^{(2)}(h,(a_1,a_2)) = \frac{2 \pi \text{erf}\left(\frac{\text{a2} h}{\sqrt{2}}\right) T(h,\text{a1})+\arctan\left(\frac{\text{a1} \text{a2}}{\sqrt{\text{a1}^2+\text{a2}^2+1}}\right) \text{erfc}\left(\frac{h \sqrt{\text{a1}^2+\text{a2}^2+1}}{\sqrt{2}}\right)}{4 \pi } + \frac{h \sqrt{1+a_2^2}}{\pi^{3/2} 2^{3/2}} \int\limits_0^{arccosh(\sqrt{\frac{1+a_1^2+a_2^2}{1+a_2^2}})} \sinh(\theta) \cdot \arctan\left(a_2 \frac{\sinh(\theta)}{\cosh(\theta)} \right) \cdot e^{-\frac{h^2}{2} (1+a_2^2) \cosh(\theta)^2} d\theta \end{eqnarray}
In particular for $h=0$ we have : \begin{equation} T^{(2)}(0,(a_1,a_2)) = \frac{1}{4 \pi} \arctan\left(\frac{\text{a1} \text{a2}}{\sqrt{\text{a1}^2+\text{a2}^2+1}}\right) \end{equation} as it should be (see An integral involving error functions and a Gaussian).