"Find a representation [using Mellin transform] of 𝑓(𝜇,𝛽) as Gauss hypergeometric function in variable 𝜇"

This is a follow-up to the first comment (by Nemo) to the posting Compute the two-fold partial integral, where the three-fold full integral is known . (I also just asked this as a comment to that question itself--but suspected that it would remain rather obscure in that form.)

After a simplifying reparameterization (transforming $$\beta$$ into $$\frac{b-1}{3}$$) suggested by Matt F., the originally-posed problem took the form, $$$$\int_{p=0}^1 \int_{q=0}^{1-p} (\mu p q (1-p-q))^b (\mu^2 q +p)^{-b-1} dq dp,$$$$ which has been employed since.

Answers to that question have, in fact, been posted by Martin Rubey and by me, but the possible one that was apparently raised by Nemo--in terms of the original $$\beta, p_{11,22}$$ (rather than $$b,p,q$$) parameterization--has not so far been fully detailed. (In his comment, Nemo wrote "Did you try calculating Mellin transform of π(π,π½) wrt to π (the double integral over $$p_{11,22}$$ can be calculated using Dirichlet's Beta integral) and then recover π(π,π½) by inverse Mellin transform (which reduces to Barnes integral in this case? I did some calculations and if not mistaken this will lead to representation of π(π,π½) as Gauss hypergeometric function in variable π".)

Nemo also later commented "MartinRubey--my observation above about this function having a Gauss hypergeometric form is consistent with this logarithmic terms because 2 parameters of this hypergeometric function coincide", but this has not yet been expanded into an answer after a request of Rubey to do so.

So, to reiterate, I would like to explicitly know the presumed representation of π(π,π½) as Gauss hypergeometric function in variable π, as this might be helpful in the pursuit of the program to construct "separability functions" presented in https://arxiv.org/abs/1701.01973.

In fact, we have been able to find, as already detailed in the earlier answers, a Gauss hypergeometric function expression for the linear (in $$\log{\mu}$$) term, $$$$w(b,\mu)=\frac{\sqrt{\pi } 4^{-b} \mu^b \left(\mu^2-1\right)^{-2 b-1} \Gamma (b+1) \, _2F_1\left(-b,-b;1;\mu^2\right)}{\Gamma \left(b+\frac{3}{2}\right)},$$$$ but not yet for the constant term'' $$v(b,\mu)$$ of the complete functional expression $$$$v(b,\mu) + w(b,\mu) \log(\mu)=$$$$ $$$$\frac{1}{\Gamma \left(b+\frac{3}{2}\right)} \sqrt{\pi } 4^{-b} \mu^b \left(\mu^2-1\right)^{-2 b-1} \Gamma (b+1) \left(\log (\mu) \sum _{k=0}^b \mu ^{2 k} \binom{b}{k}^2-\left(\mu^2-1\right) \sum _{k=1}^b \mu^{2 k-2} \sum _{i=0}^{k-1} \binom{b}{i}^2 (\psi ^{(0)}(b-i+1)-\psi ^{(0)}(i+1))\right),$$$$ where, $$$$\sum _{k=0}^b \mu ^{2 k} \binom{b}{k}^2=\, _2F_1\left(-b,-b;1;\mu ^2\right).$$$$

In trying to implement the strategy of Nemo, I performed the integration of $$$$\int_{p=0}^1 \int_{q=0}^{1-p} (\mu p q (1-p-q))^b (\mu^2 q +p)^{-b-1} dq dp,$$$$ for $$b =1,\ldots,5$$, then requested Mathematica to take the Mellin transform with respect to $$\mu$$ of the results.
For $$b=1$$, I obtained $$$$\frac{1}{3} \left(\mathcal{M}_\mu\left[\frac{\mu}{\left(\mu^2-1\right)^3}\right](s)+\mathcal{M}_\mu\left[\frac{\mu \log (\mu)}{\left(\mu^2-1\right)^3}\right](s)-\mathcal{M}_\mu\left[\frac{\mu^3}{\left(\mu^2-1\right)^ 3}\right](s)+\mathcal{M}_\mu\left[\frac{\mu^3 \log (\mu)}{\left(\mu^2-1\right)^3}\right](s)\right),$$$$ for $$b=2$$, $$$$\frac{1}{30} \left(3 \mathcal{M}_\mu\left[\frac{\mu^2}{\left(\mu^2-1\right)^5}\right](s)-3 \mathcal{M}_\mu\left[\frac{\mu^6}{\left(\mu^2-1\right)^5}\right](s)-\frac{\pi (s-2) s \left(e^{i \pi s} \left(\pi s^2-2 (\pi -2 i) s-4 i\right)-4 i (s-1)\right)}{32 \left(-1+e^{i \pi s}\right)^2}\right),$$$$ for \$b=3, $$\frac{\pi (s-3) (s-1) (s+1) \left(2 i (3 (s-2) s-1)+e^{i \pi s} (6 i (s-2) s+\pi (s-3) (s-1) (s+1)-2 i)\right)}{161280 \left(1+e^{i \pi s}\right)^2} +$$ $$-\frac {11} {420} \mathcal {M} _\mu [\frac {\mu^9} { (\mu^2 - 1 )^7} ] (s) - \frac {9} {140} \mathcal {M} _\mu [\frac {\mu^7} { (\mu^2 - 1 )^7} ] (s) + \frac {9} {140} \mathcal {M} _\mu [\frac {\mu^5} { (\mu^2 - 1 )^7} ] (s) + \frac {11} {420} \mathcal {M} _\mu [\frac {\mu^3} { (\mu^2 - 1 )^7} ] (s),$$ ....
So, since the required Mellin transforms do not all seem explicitly implementable, it appears to me that the Nemo goal of obtaining a (complete) representation of $$f(\mu,b)$$ (using the Matt. F. parameters) does not seem realizable.
The integral can be solved using the Mellin representation of a sum, which is $$(A+B)^{-s} = \int_C \frac{\Gamma(s-t) \Gamma(t)}{\Gamma(s)} A^{-t} B^{t-s} \frac{\mathrm{d} t}{2 \pi \mathrm{i}} \,,$$ valid (at least) for $$\Re A, \Re B, \Re s > 0$$, and where the integration contour $$C$$ is the usual one.
We write $$(\mu^2 q + p)^{-(b+1)} = \int_C \frac{\Gamma(b+1-t) \Gamma(t)}{\Gamma(b+1)} (\mu^2 q)^{-t} p^{t-(b+1)} \frac{\mathrm{d} t}{2 \pi \mathrm{i}} \,,$$ and then interchange integrations (justified in the usual way) to obtain $$\int_0^1 \int_0^{1-p} (\mu p q (1-p-q))^b (\mu^2 q + p)^{-b-1} \mathrm{d}q \, \mathrm{d}p = \int_C \mu^{b-2t} \frac{\Gamma(b+1-t) \Gamma(t)}{\Gamma(b+1)} \int_0^1 \int_0^{1-p} (1-p-q)^b q^{b-t} p^{t-1} \mathrm{d}q \, \mathrm{d}p \, \frac{\mathrm{d} t}{2 \pi \mathrm{i}} \,.$$ The integral over $$q$$ gives $$\int_0^{1-p} (1-p-q)^b q^{b-t} \mathrm{d}q = (1-p)^{2b+1-t} \frac{\Gamma(1+b) \Gamma(1+b-t)}{\Gamma(2+2b-t)} \,,$$ and the remaining integral over $$p$$ results in $$\int_0^1 p^{t-1} (1-p)^{2b+1-t} \mathrm{d}p = \frac{\Gamma(2+2b-t) \Gamma(t)}{\Gamma(2+2b)} \,,$$ such that $$\int_0^1 \int_0^{1-p} (\mu p q (1-p-q))^b (\mu^2 q + p)^{-b-1} \mathrm{d}q \, \mathrm{d}p = \int_C \mu^{b-2t} \frac{\Gamma^2(b+1-t) \Gamma^2(t)}{\Gamma(2+2b)} \frac{\mathrm{d} t}{2 \pi \mathrm{i}} \,.$$ The remaining inverse Mellin transform finally gives a hypergeometric function, and we obtain the result $$\int_0^1 \int_0^{1-p} (\mu p q (1-p-q))^b (\mu^2 q + p)^{-b-1} \mathrm{d}q \, \mathrm{d}p = \frac{\pi \mu^b \, \Gamma^2(b+1)}{4^{1+2b} \, \Gamma^2\left( \frac{3}{2} + b \right)} {}_2\mathrm{F}{}_1\left( b+1, b+1; 2 (b+1); 1-\mu^2 \right) \,.$$ The confluent hypergeometric identity https://functions.wolfram.com/HypergeometricFunctions/Hypergeometric2F1/17/02/08/0003/ with $$n = m = 0$$ finally turns this into $$\int_0^1 \int_0^{1-p} (\mu p q (1-p-q))^b (\mu^2 q + p)^{-b-1} \mathrm{d}q \, \mathrm{d}p = - \frac{\pi \mu^b \, \Gamma(2b+2)}{2^{1+4b} \, \Gamma^2\left( \frac{3}{2} + b \right)} \left[ {}_2\mathrm{F}{}_1\left( b+1, b+1; 1; \mu^2 \right) \ln \mu - \sum_{k=0}^\infty \frac{\Gamma^2(b+1+k)}{\Gamma^2(b+1) \Gamma^2(k+1)} \mu^{2k} \left( \psi(k+1) - \psi(b+k+1) \right) \right] \,,$$ which shows explicitly the logarithmic terms. For each integer $$b$$, the hypergeometric function becomes a rational function of $$\mu^2$$, and the summand becomes a polynomial in $$k$$, which after summation gives another rational function of $$\mu^2$$. I don't think that there is a representation of the sum using GauΓ hypergeometric functions, but possibly of generalized ones.