Compute the two-fold partial integral, where the three-fold full integral is known I have the following trivariate ($\rho_{11}, \rho_{22}, \mu$) function
\begin{equation}
4 \mu ^{3 \beta +1} \rho_{11}^{3 \beta +1} \left(-\rho_{11}-\rho_{22}+1\right){}^{3 \beta +1}
   \rho_{22}^{3 \beta +1} \left(\mu ^2 \rho_{22}+\rho_{11}\right){}^{-3 \beta -2},
\end{equation}
the (three-fold) integral (for $\beta$ nonnegative integer) of which over $\mu \in [0,1]$, $\rho_{11} \in [0,1]$,$\rho_{22} \in [0,1- \rho_{11}]$ is 
\begin{equation}
\frac{\Gamma \left(\frac{3 \beta }{2}+1\right)^4}{\Gamma (6 \beta +4)}.
\end{equation}
I would like to know the (two-fold) integral, say $f(\mu,\beta)$, over $\rho_{11} \in [0,1]$,$\rho_{22} \in [0,1- \rho_{11}]$. (The three-fold integral result is based on the application of the Mathematica FindSequenceFunction command, and does not constitute a formal proof--as remarks of C. Dunkl lead me to state.)
A reference for the (quantum-information-theoretic) background of this problem is sec. II.A.1 of Slater - Extended Studies of Separability Functions and Probabilities and the Relevance of Dyson Indices, in particular eq. (8) there. (For later related work, see secs. III,IV of Slater - Master Lovas-Andai and Equivalent Formulas Verifying the 8/33 Two-Qubit Hilbert-Schmidt Separability Probability and Companion Rational-Valued Conjectures, where I try to relate  Lovas-Andai and Slater "separability functions".)
Mathematica computes the integral for specific nonnegative integer $\beta$, but apparently not for general $\beta$. (I can compute for $\beta=1,2,\ldots$ and then try to employ the Mathematica command FindSequenceFunction to obtain the general [two-fold] rule--which is, I think, how I got the three-fold result.)
For example, for $\beta=1$, we have
\begin{equation}
\frac{\mu ^4 \left(12 \left(\mu^8+16\mu ^6+36 \mu ^4+16\mu
   ^2+1\right) \log (\mu )-5 \left(5 \mu ^8+32 \mu ^6-32 \mu ^2-5\right)\right)}{945
   \left(\mu ^2-1\right)^9}
\end{equation}
and for $\beta=2$,
\begin{equation}
\frac{\mu ^7
\left(140 \left(\mu^2+1\right)\left(\mu ^{12}+48 \mu ^{10}+393 \mu ^8+832 \mu ^6+393 \mu ^4+48 \mu ^2+1\right)\log(\mu )
-(\mu^2-1) \left(363 \mu ^{12}+10310 \mu ^{10}+58673 \mu^8+101548 \mu ^6+58673 \mu ^4+10310 \mu ^2+363\right)\right)
}{900900 \left(\mu ^2-1\right)^{15}}.
\end{equation}
 A: Using FriCAS it is actually not hard to guess the complete solution.  However, this is not a proof, one would have to show


*

*that the integral is indeed a linear polynomial of $\log\mu$, (@Nemo?)

*that both the linear and the constant term of this polynomial must satisfy a second order recurrence with polynomial coefficient (@Nemo?),

*and that the degrees of the coefficients in the recurrence can be bounded.


It came as a surprise, but it probably shouldn't, that the recurrences for both the linear and he constant term are the same, they differ only in the initial conditions.
Anyway, here is code and solution.  Be warned that I took out the $\mu$, and a factor of $4$ in the linear term.
    
f := (p*q*(1-p-q))^b*(m^2*q+p)^(-b-1)
)se fu ca all fint
fint(b:NNI):EXPR INT == integrate(integrate(eval(f, 'b=b), q=0..1-p, "noPole")::EXPR INT, p=0..1, "noPole")::EXPR INT
logm := first kernels fint(1)::EXPR INT
fintUP b == eval(fint(b), logm = lm)::UP(lm, FRAC POLY INT)

r1 := guessPRec([coefficient(fintUP b, 1) for b in 0..15], indexName=='b)
r2 := guessPRec([coefficient(fintUP b, 0) for b in 0..15], indexName=='b)

)expose RECOP
fun := operator 'f
r1DMP := eval(getEq(r1.1), [fun b = 'fb0, fun(b+1) = 'fb1, fun(b+2) = 'fb2])::DMP(['fb,'fb1,'fb2], POLY INT);
r2DMP := eval(getEq(r2.1), [fun b = 'fb0, fun(b+1) = 'fb1, fun(b+2) = 'fb2])::DMP(['fb,'fb1,'fb2], POLY INT);

factor coefficient(r1DMP, 'fb0, 1)
factor coefficient(r1DMP, 'fb1, 1)
factor coefficient(r1DMP, 'fb2, 1)

The recurrence for both the constant and the linear term is:
$$
-(b+1)^2 f(b) + 2 (2b+3)^2(\mu^2 +1) f(b+1) - 4 (2b+3)(2b+5)(\mu^2-1)^2 f(b+2) = 0
$$
The initial conditions for the linear term are:
$$
f(0)=\frac{1}{2(\mu^2-1)},\quad f(1)=\frac{\mu^2+1}{12(\mu^2-1)^3}
$$
The initial conditions for the constant term are:
$$
f(0)=0,\quad f(1)=\frac{1}{3(\mu^2-1)^2}
$$
A: The most concise answer to the problem (in the Matt F. streamlined form)
\begin{equation}
 \int_{p=0}^1 \int_{q=0}^{1-p} (\mu p q (1-p-q))^b (\mu^2 q +p)^{-b-1} dq dp,
\end{equation} 
so far given is
\begin{equation}
\frac{\pi  (-1)^{2 b} 4^{-2 b-1} \mu^{-b-2} \Gamma (b+1)^2 \, _2F_1\left(b+1,b+1;2
   (b+1);1-\frac{1}{\mu^2}\right)}{\Gamma \left(b+\frac{3}{2}\right)^2}.
\end{equation}
Its validity can be seen by, first, performing the inner integration of $q \in [0,1-p]$, obtaining thereby,
\begin{equation}
\int_{p=0}^1-\frac{(p-1)^{2 b+1} \mu^b \Gamma (b+1)^2 \, _2F_1\left(b+1,b+1;2 (b+1);\frac{(p-1)
   \mu^2}{p}\right)}{p \Gamma (2 (b+1))} dp.
\end{equation}
Then, the hypergeometric argument $\frac{(p-1)
   \mu^2}{p}$ is transformed to $v$, by performing the change-of-variables $p\to \frac{\mu^2}{\mu^2-v}$. The one-fold integration problem, which can be carried out, then, becomes
\begin{equation}
\int_{v=-\infty}^0 \mu ^b \left(-v^{2 b+1}\right) \Gamma (b+1)^2 \, _2\tilde{F}_1(b+1,b+1;2 b+2;v) \left(\mu
   ^2-v\right)^{-2 (b+1)} dv,
\end{equation}
yielding us the indicated answer given above here to the two-fold integration (in the Matt F. form).
However, to perform the third step of the three-fold integration, that is,
\begin{equation}
\int_{\mu=0}^1 \frac{\pi  (-1)^{2 b} 4^{-2 b-1} \mu^{-b-2} \Gamma (b+1)^2 \, _2F_1\left(b+1,b+1;2
   (b+1);1-\frac{1}{\mu^2}\right)}{\Gamma \left(b+\frac{3}{2}\right)^2} d \mu,
\end{equation}
we first perform--as pointed out by Charles Dunkl--the quadratic transformation 15.8.13 given in 
https://dlmf.nist.gov/15.8#iii   . That is 
\begin{equation}
\, _2F_1(a,b;2 b;z) = (1-\frac{1}{2} z)^{-a} \, _2F_1(\frac{1}{2} a,\frac{1}{2} a+\frac{1}{2};b+\frac{1}{2} ;(\frac{z}{2-z})^2)
\end{equation}
on the integrand and make the change-of-variables $\mu\to \sqrt{\frac{2}{\sqrt{Y}+1}-1}$ and now integrate over $Y \in [1,0]$, that is,
\begin{equation}
\int_{Y=1}^0 \frac{\pi  2^{-4 b-3} e^{2 i \pi  b}
   \left(\frac{2}{\sqrt{Y}+1}-1\right)^{\frac{1-b}{2}}
   \left(\frac{1}{1-\sqrt{Y}}\right)^{-b} \Gamma (b+1)^2 \,
   _2F_1\left(\frac{b+1}{2},\frac{b+2}{2};b+\frac{3}{2};Y\right)}{\left(Y-\sqrt{Y}\right) \Gamma \left(b+\frac{3}{2}\right)^2} dY.
\end{equation} 
The result of this integration is 
\begin{equation}
\frac{e^{2 i \pi  b} \Gamma \left(\frac{b+1}{2}\right)^4}{4 \Gamma (2 b+2)}.
\end{equation}
For integral $b$, the phase factor $e^{2 i \pi  b}$ is, of course, 1, and for half-integral $b$, it is -1.
Transforming back to the original (Dyson-index-like parameter) $\beta$, using $b \to 1+3 \beta$, we obtain
\begin{equation}
\frac{e^{6 i \pi  \beta } \Gamma \left(\frac{3 \beta }{2}+1\right)^4}{4 \Gamma (6
   \beta +4)}.
\end{equation}
This agrees (up to the omitted--in the Matt F. parameterization--factor 4) with the original assertion of 
\begin{equation}
\frac{\Gamma \left(\frac{3 \beta }{2}+1\right)^4}{\Gamma (6 \beta +4)},
\end{equation}
for integral  $\beta$.
Here is a link--through the wolframcloud--to a Mathematica notebook presenting the main later stages of our argument.
https://www.wolframcloud.com/objects/e8bd5d73-06a4-4798-8485-8e2597f660eb
