I want to compute the integral of the product of two inverse regularized incomplete beta functions over $[0,1]$ in closed form, that is, to evaluate $$ J = \int_0^1 I_t^{-1}(a_1,b_1) \: I_t^{-1}(a_2,b_2) \: \mathrm{d}t $$ in terms of parameters $a_1, b_1, a_2, b_2 > 0$.
Let us first collect few useful facts:
(Indefinite integral) $\displaystyle\int I_{t}^{-1}(a,b)\:\mathrm{d}t = \frac{\left(I_t^{-1}(a,b)\right)^{a+1}}{(a+1)\:B(a,b)}\:_2 F_1\left(a+1, 1-b; a+2; I_t^{-1}(a,b)\right) + \text{constant}$.
(Derivative) $\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}I_t^{-1}(a,b) = B(a,b)\:\left(I_t^{-1}(a,b)\right)^{1-a}\:\left(1-I_t^{-1}(a,b)\right)^{1-b}$.
(Boundary evaluation) $I_0^{-1}(a,b) = 0, \quad I_1^{-1}(a,b)=1$.
(Gauss's Hypergeometric theorem) $_2 F_{1}(p,q;r;1) = \dfrac{\Gamma(r) \:\Gamma(r-p-q)}{\Gamma(r-p)\:\Gamma(r-q)}$.
(Complete Beta and Gamma functions) $B(a,b) = \dfrac{\Gamma(a)\Gamma(b)}{\Gamma(a+b)}$.
Integrating-by-parts the original integral $J$ by taking $I_t^{-1}(a_1,b_1)$ as the first function, and $I_t^{-1}(a_2,b_2)$ as the second, we are led to $J = J_1 - J_2$, with a simple expression for \begin{align} J_1 & = \left. \left[I_t^{-1}(a_1,b_1) \int I_t^{-1}(a_2,b_2) \: \mathrm{d}t \right] \right|_{t=0}^{t=1} \\[10pt] & = \frac{a_2}{a_2 + b_2}, \quad \text{(using facts 1,3,4,5, and } \Gamma(z+1) = z\Gamma(z)), \end{align} and a complicated-looking expression for \begin{align} J_2 & = \int_0^1 \left(\frac{\mathrm{d}}{\mathrm{d}t}I_t^{-1}(a_1,b_1) \int I_t^{-1}(a_2,b_2)\:\mathrm{d}t\right)\:\mathrm{d}t \\[10pt] & = \frac{B(a_1,b_1)}{(a_2+1) B(a_2,b_2)} \int_0^1 \left(I_t^{-1}(a_1, b_1)\right)^{1-a_1} \left(1 - I_t^{-1}(a_1,b_1)\right)^{1-b_1} \left(I_t^{-1}(a_2,b_2)\right)^{a_2+1}\:_2 F_1\left(a_2+1,1-b_2;a_2+2;I_t^{-1}(a_2,b_2)\right)\mathrm{d}t, \quad\text{(using facts 1,2)} \end{align}
Can $J_2$ be further simplified? Or perhaps there is a more elegant way to compute $J$ than my naive attempt via integration-by-parts? Any help, suggestions or ideas are welcome.
