# 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 $$$$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},$$$$ 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 $$$$\frac{\Gamma \left(\frac{3 \beta }{2}+1\right)^4}{\Gamma (6 \beta +4)}.$$$$ 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}]$$.

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 $$$$\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}$$$$ and for $$\beta=2$$, $$$$\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}}.$$$$

• Did you try calculating Mellin transform of $f(\mu,\beta)$ wrt to $\mu$ (the double integral over $p_{11,22}$ can be calculated using Dirichlet's Beta integral) and then recover $f(\mu,\beta)$ 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 $f(\mu,\beta)$ as Gauss hypergeometric function in variable $\mu$. – Nemo Feb 11 at 15:24
• I would ask for $$\int_{p=0}^1\int_{q=0}^{1-p}(\mu pq(1-p-q))^b(\mu^2q+p)^{-b-1}dp\ dq.$$ The messy variable names from the original are valuable context, but can be left out of the integral for clarity. – Matt F. Feb 11 at 15:57
• It seems that the @MattF.'s integral is a linear polynomial in $\log\mu$ (I did not attempt to prove this), whose constant term has denominator $(\mu^2-1)^{2b}$ and whose linear term has denominator $(\mu^2-1)^{2b+1}$. – Martin Rubey Feb 11 at 18:28
• The linear polynomials (in log $\mu$) that Martin Rubey observes appear to be multiples of $\mu ^{3 \beta +1} \left(\mu ^2-1\right)^{-6 \beta -3}$. – Paul B. Slater Feb 11 at 19:28
• In fact, the linear term (in @MattF.'s integral) is apparently $\frac{1}{(m^2-1)^{2b+1}}\frac{1}{4(4b^2-1)}\frac{b}{\binom{2(b-1)}{b-1}}\sum_k\binom{b}{k}^2m^{2k}$. – Martin Rubey Feb 11 at 19:45

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}$$

• That's great, Martin! I'll have to cogitate on this further in the A. M. First I heard of FriCAS. Can we get an explicit expression for the constant term akin to that I have for the linear (in log(\mu)) term? (I was getting pretty confident that the two terms had "a lot in common".) Also, you've been apparently using $m$ and $\mu$ for the same variable. Once things are clearly settled, I'll have to think about how to employ the results in the quantum-information-theoretic context that I referenced in my initial question. – Paul B. Slater Feb 14 at 4:59
• So, Martin, I presume there is a summation formula for the constant term, parallel to that you gave for the linear term. Mathematica "automatically" converted that formula into the hypergeometric-based one. So, a parallel formula should give me the natural companion v(b,\mu) to my w(b,\mu) expression. – Paul B. Slater Feb 14 at 13:13

To begin, a partial/half answer based on the summation formula of Martin Rubey, given that the two-fold integration result in the Matt F. streamlined reformulation takes the form $$$$v(b,\mu) + w(b,\mu) \log(\mu) = \int_{p=0}^1 \int_{q=0}^{1-p} (\mu p q (1-p-q))^b (\mu^2 q +p)^{-b-1} dq dp,$$$$ is that (in line with the "hypergeometric" comment of Nemo above) $$$$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)}.$$$$ The first term in the commented Rubey summation formula, $$\frac{b \left(\mu ^2-1\right)^{-2 b-1} \, _2F_1\left(-b,-b;1;\mu ^2\right)}{4 \left(4 b^2-1\right) \binom{2 (b-1)}{b-1}}$$, for $$w(b,\mu)$$ above (with $$m$$ corrected to $$\mu$$) required an additional factor of $$4 \mu^b$$ (as Ruhey noted in a subsequent comment). (Mathematica insists, it seems, on converting his summation to hypergeometric form. It would be of interest to know how Ruhey arrived at his formula--and if it might be of potential use in obtaining $$v(b,\mu)$$, as well.)

We have, additionally, trying to address the other "half" of the problem, found that $$$$v(b,1)=\frac{\pi 4^{-2 b-1} \Gamma (b+1)^2}{\Gamma \left(b+\frac{3}{2}\right)^2}.$$$$

Further, for $$v(b,\frac{1}{2})$$, Mathematica has given DifferenceRoot[Function[{[FormalY],[FormalN]},{(1+[FormalN])^2 [FormalY][[FormalN]]-5 (3+2 [FormalN])^2 [FormalY][1+[FormalN]]+(135+144 [FormalN]+36 [FormalN]^2) [FormalY][2+[FormalN]]==0,[FormalY][1]==-(8/27),[FormalY][2]==-(8/81)}]][b]

The TeX version of this that Mathematica gives is not accepted here because the "Argument to unicode must be a number". (This appears to be a problem with the TeXForm command of Mathematica https://mathematica.stackexchange.com/questions/191414/nonacceptance-by-stackexchange-site-of-mathematica-texform-employing-unicode .)

Along similar lines, for $$v(b,1/3)$$, we have DifferenceRoot[ Function[{[FormalY], [FormalN]}, {9 (1 + [FormalN])^2 \ [FormalY][[FormalN]] - 60 (3 + 2 [FormalN])^2 [FormalY][ 1 + [FormalN]] + (3840 + 4096 [FormalN] + 1024 [FormalN]^2) [FormalY][2 + [FormalN]] == 0, [FormalY][1] == -(9/64), [FormalY][2] == -(81/4096)}]][b]

So, "DifferenceRoot" formulas appear to be the case for other specific values of $$v(b,\mu)$$ with $$\mu \neq 1$$.

Perhaps, we will attempt ("manually") to present the DifferenceRoot formulas above in more standard, proper TeX form (per comment of Somos in https://mathematica.stackexchange.com/questions/191414/nonacceptance-by-stackexchange-site-of-mathematica-texform-employing-unicode ).

It may be helpful to note that the inner integral of the Matt F. reformulation can be performed, yielding $$$$\int_{q=0}^{1-p} (\mu p q (1-p-q))^b (\mu^2 q +p)^{-b-1} dq =-\frac{(p-1) \Gamma (b+1)^2 \left((p-1)^2 \mu\right)^b \, _2\tilde{F}_1\left(b+1,b+1;2 (b+1);\frac{(p-1) \mu^2}{p}\right)}{p},$$$$ where the regularized hypergeometric function is indicated.

So, we would like the counterpart for $$v(b,\mu)$$ of the Rubey formula for $$w(b,\mu)$$, that is, $$$$4 u^b \left(u^2-1\right)^{-2 b-1} \frac{1}{4 \left(4 b^2-1\right)} \frac{b} {\binom{2 (b-1)}{b-1}} \Sigma_{k=0}^b u^{2 k} \binom{b}{k}^2,$$$$ which he apparently obtained using the general purpose computer algebra system, FriCAS.

A new, interesting observation is that if we multiply $$v(b,\mu)$$ by $$$$(-1+\mu^2)^{2 b+1} \mu^{-b} P,$$$$ where, with $$\gamma$$ being Euler's constant, $$\psi ^{(0)}(b+1)$$ being the zero-th derivative of the digamma function, and $$H_b$$, the $$b$$-th harmonic number, $$$$P= \frac{\sqrt{\pi } 2^{-2 b} b! (\psi ^{(0)}(b+1)+\gamma )}{\left(b+\frac{1}{2}\right)!} = \frac{\sqrt{\pi } 4^{-b} b! H_b}{\Gamma \left(b+\frac{3}{2}\right)},$$$$ we obtain ("asymmetric"-type) polynomials in $$\mu$$ of degree $$2 b$$, all with constant terms equal to 1 and highest terms equal to $$-u^{2 b}$$. For example, for $$b=1$$, we have $$$$1-\mu^2$$$$ for $$b=2$$, $$$$1-\mu^4$$$$ for $$b=3$$, $$$$-\mu ^6-\frac{27 \mu ^4}{11}+\frac{27 \mu ^2}{11}+1$$$$ for $$b=4$$, $$$$-\mu ^8-\frac{32 \mu ^6}{5}+\frac{32 \mu ^2}{5}+1,$$$$ for $$b=5$$, $$$$-\mu ^{10}-\frac{1625 \mu ^8}{137}-\frac{2000 \mu ^6}{137}+\frac{2000 \mu ^4}{137}+\frac{1625 \mu ^2}{137}+1, \ldots$$$$ So, one needs to find the governing rule for these polynomials, which would complete the formula for $$v(b,\mu)$$.

• I am quite sure that I do get the factor $(\mu^2-1)^{2b+1}$, could you please recheck? For example, for $b=2$, I obtain $$\frac{\mu^4+4\mu^2+1}{15(\mu^2-1)^5}$$. – Martin Rubey Feb 13 at 6:59
• Thanks, Martin! I did some quick checking this morning--before off to campus--and I can't, at the moment, get anything to agree. Seems that I may need a factor of u^b/(u^2-1)^(2 b+1) in my formula. Embarrassing! – Paul B. Slater Feb 13 at 14:44
• Oh yes, I took out the $\mu^b$ from the integral! – Martin Rubey Feb 13 at 15:29
• So, the revised answer accounts for these last matters. – Paul B. Slater Feb 13 at 22:37
• Matt F. I think the result (according to Mathematica) should be $\, _2F_1\left(-b,-b;1;u^2\right)$. $k$ is the index of summation and should not appear in the answer. – Paul B. Slater Feb 14 at 21:48