# Evaluate over a two-dimensional domain, the integral of (hypergeometric-based) $f(d,k)$, the result for $f(d,0)$ being known

I view this as both a mathematics and a Mathematica question -- so apologies if it is thought I should have sent it alternatively to the mathematics stack exchange.

I want to perform the two-dimensional integration (or, possibly, reduce to a one-dimensional integration)

Integrate[Y^(-1 + d)
Hypergeometric2F1Regularized[d/2, -k, (2 + d)/
2, ((Y^2 - e^2 Subscript[r, 14]^2) (-1 + e^2 Subscript[r, 14]^2))/(
e^2 (Y^2 - Subscript[r, 14]^2) (-1 + Subscript[r, 14]^2))] (1/(
e Subscript[r, 14]))^(
1 + d) (1 + Y^2 (1 - 1/Subscript[r, 14]^2) - Subscript[r, 14]^2)^
k (1 - e^2 Subscript[r, 14]^2)^(
d/2) (-Y^2 + e^2 Subscript[r, 14]^2)^(d/2), {Subscript[r, 14], 0,
1}, {Y, e Subscript[r, 14]^2, e Subscript[r, 14]},
Assumptions -> d >= 1 && k >= 0 && 0 < e <= 1]


So, $$d$$ and $$k$$ are parameters, and $$Y$$ and $$r_{14}$$ the variables of integration, with $$e$$ being a free variable.

In $$\mathrm\TeX$$, the integrand is the product of

$$\begin{equation} Y^{d-1} \left(\frac{1}{r_{14} \epsilon }\right){}^{d+1} \left(1-r_{14}^2 \epsilon ^2\right){}^{d/2} \left(\left(1-\frac{1}{r_{14}^2}\right) Y^2-r_{14}^2+1\right){}^k \left(r_{14}^2 \epsilon ^2-Y^2\right){}^{d/2} \end{equation}$$

and

$$\begin{equation} \, _2\tilde{F}_1\left(\frac{d}{2},-k;\frac{d+2}{2};\frac{\left(r_{14}^2 \epsilon ^2-1\right) \left(Y^2-r_{14}^2 \epsilon ^2\right)}{\left(r_{14}^2-1\right) \epsilon ^2 \left(Y^2-r_{14}^2\right)}\right) . \end{equation}$$

The two-dimensional domain of integration is

$$\begin{equation} r_{14} \in [0,1], \hspace{.25in} Y \in [\varepsilon r_{14}, \varepsilon^2 r_{14}] . \end{equation}$$

For $$k=0$$, the integral evaluates (as can be confirmed by setting $$d$$ to a positive integer -- even integers evaluate more readily) to

1/4 e^(-1 + d)
Gamma[d/2] Gamma[
d] HypergeometricPFQRegularized[{-(d/2), d/2, d}, {1 + d/2,
1 + (3 d)/2}, e^2]


That is, $$\begin{equation} \frac{1}{4} \epsilon ^{d-1} \Gamma \left(\frac{d}{2}\right) \Gamma (d) \, _3\tilde{F}_2\left(-\frac{d}{2},\frac{d}{2},d;\frac{d}{2}+1,\frac{3 d}{2}+1;\epsilon ^2\right). \end{equation}$$

I've been trying all possible integration-by-parts combinations with no success to this point in time.

For even, nonnegative $$d$$ and nonnegative $$k$$, the integrand evaluates to a polynomial in $$e$$. For odd $$d$$, logs and polylogs appear.

The question stated here pertains to the issue discussed in sec. IX.B of my posting, Qubit-qudit separability/PPT-probability investigations, including Lovas-Andai formula advancements, of finding an "extended Lovas-Andai master formula", denoted there by $$\tilde{\chi}_{d,k}(\varepsilon)$$.

Charles F. Dunkl has provided me with essentially the answer to the bivariate ($r_{14}, r_{23}$) integration question at hand. (However, I had brought to his attention, a problem with a slightly reexpressed integrand, using a change-of-variables, now designated $x,y$. Also, I had omitted certain accompanying [gamma,...] factors, not functions of the two variables, from the stated problem for the sake of conciseness. These are now present in his write-up. The important parameters $d,k$ and (not integrated out) variable $\varepsilon$ remain. (Let me note that CFD works in Maple.)
I now provide his write-up (certain of the TeX commands do not seem to "take" in this stack-exchange setting). The central problem is the evaluation of $I(\varepsilon)$, which is given near the end.
Let us simplify \begin{align*} I\left( \varepsilon\right) & :=\frac{\Gamma\left( 1+d+k\right) ^{2}% }{\Gamma\left( \frac{d}{2}\right) ^{3}\Gamma\left( 1+\frac{d}{2}+k\right) \Gamma\left( 1+k\right) }\times\frac{2}{d}\varepsilon^{-d}\\ & \times\int_{0}^{1}dx\int_{\varepsilon^{2}x}^{\varepsilon^{2}}dy[\left\{ \left( 1-x\right) \left( 1-y\right) \right\} ^{k}\left( xy\right) ^{d/2-1}\left\{ \left( \varepsilon^{2}-y\right) \left( 1-x\varepsilon ^{2}\right) \right\} ^{d/2}\\ & \times~_{2}F_{1}\left( -k,\frac{d}{2};1+\frac{d}{2};T\right) ]\\ T & :=\frac{\left( \varepsilon^{2}-y\right) \left( 1-x\varepsilon ^{2}\right) }{\left( 1-x\right) \left( 1-y\right) \varepsilon^{2}}. \end{align*} First we apply the transformation $_{2}F_{1}\left( a,b;c;t\right) =\left( 1-t\right) ^{-a}~_{2}F_{1}(a,c-b;c;\frac{t}{t-1}$), but the series on the right side only converges if $a$ is a negative integer or $t<\frac{1}{2}$, not the case in our application, thus \textbf{henceforth assume} $k=0,1,2,3,\ldots$then% \begin{align*} 1-T & =\frac{\left( 1-\varepsilon^{2}\right) \left( y-x\varepsilon ^{2}\right) }{\left( 1-x\right) \left( 1-y\right) \varepsilon^{2}},\\ \frac{T}{T-1} & =-\frac{\left( \varepsilon^{2}-y\right) \left( 1-x\varepsilon^{2}\right) }{\left( 1-\varepsilon^{2}\right) \left( y-x\varepsilon^{2}\right) }; \end{align*} the integrand becomes% \begin{align*} & \left( xy\right) ^{d/2-1}\left( 1-\varepsilon^{2}\right) ^{k}% \varepsilon^{-2k}\left( y-x\varepsilon^{2}\right) ^{k}\left( 1-x\varepsilon ^{2}\right) ^{d/2}\left( \varepsilon^{2}-y\right) _{~}^{d/2}\\ & \times~_{2}F_{1}\left( -k,1;1+\frac{d}{2};-\frac{\left( \varepsilon ^{2}-y\right) \left( 1-x\varepsilon^{2}\right) }{\left( 1-\varepsilon ^{2}\right) \left( y-x\varepsilon^{2}\right) }\right) . \end{align*} Substitute $y=\varepsilon^{2}u$ so $dy=\varepsilon^{2}du$ and $0\leq x\leq u\leq1.$ This gives a factor of $\varepsilon^{2d}$ in front of% \begin{align*} & \int\limits_{0\leq x\leq u\leq1}\int dx~du~\left( xu\right) ^{d/2-1}\\ & \sum_{j=0}^{k}\frac{\left( -k\right) _{j}}{\left( 1+\frac{d}{2}\right) _{j}}\left( -1\right) ^{j}\left( 1-x\varepsilon^{2}\right) ^{d/2+j}\left( 1-u\right) ^{d/2+j}\left( 1-\varepsilon^{2}\right) ^{k-j}\left( u-x\right) ^{k-j}. \end{align*} Isolate the $x$-integral (use the negative binomial series for $\left( 1-x\varepsilon^{2}\right) ^{d/2+j}$)% \begin{align*} & \int_{0}^{u}x^{d/2-1}\left( u-x\right) ^{k-j}\sum_{i=0}^{\infty}% \frac{\left( -\frac{d}{2}-j\right) _{i}}{i!}x^{i}\varepsilon^{2i}dx\\ & =\sum_{i=0}^{\infty}\frac{\left( -\frac{d}{2}-j\right) _{i}}{i!}% \frac{\Gamma\left( \frac{d}{2}+i\right) \Gamma\left( k-j+1\right) }% {\Gamma\left( \frac{d}{2}+i+k-j+1\right) }\varepsilon^{2i}u^{d/2+i+k-j}, \end{align*} by use of $\int_{0}^{u}x^{\alpha-1}\left( u-x\right) ^{\beta-1}% dx=u^{\alpha+\beta-1}B\left( \alpha,\beta\right)$. The inner $u$-integral is% [ \int_{0}^{1}u^{d/2-1}u^{d/2+i+k-j}\left( 1-u\right) ^{d/2+j}du=\frac {\Gamma\left( d+i+k-j\right) \Gamma\left( \frac{d}{2}+j+1\right) }% {\Gamma\left( \frac{3d}{2}+i+k+1\right) }. ] Thus the integral is% \begin{align*} & \sum_{j=0}^{k}\frac{\left( -k\right) _{j}}{\left( 1+\frac{d}{2}\right) _{j}}\varepsilon^{2d}\left( 1-\varepsilon^{2}\right) ^{k-j}\left( -1\right) ^{j}\\ & \times\sum_{i=0}^{\infty}\frac{\left( -\frac{d}{2}-j\right) _{i}}{i!}% \frac{\Gamma\left( \frac{d}{2}+i\right) \Gamma\left( k-j+1\right) }% {\Gamma\left( \frac{d}{2}+i+k-j+1\right) }\varepsilon^{2i}\frac {\Gamma\left( d+i+k-j\right) \Gamma\left( \frac{d}{2}+j+1\right) }% {\Gamma\left( \frac{3d}{2}+i+k+1\right) }\\ & =\sum_{j=0}^{k}\frac{\left( -k\right) _{j}}{\left( 1+\frac{d}{2}\right) _{j}}\left( -1\right) ^{j}\varepsilon^{2d}\left( 1-\varepsilon^{2}\right) ^{k-j}\frac{\Gamma\left( \frac{d}{2}\right) \Gamma\left( k-j+1\right) \Gamma\left( d+k-j\right) \Gamma\left( \frac{d}{2}+j+1\right) }% {\Gamma\left( \frac{d}{2}+k-j+1\right) \Gamma\left( \frac{3d}% {2}+k+1\right) }\\ & \times\sum_{i=0}^{\infty}\frac{\left( -\frac{d}{2}-j\right) _{i}\left( \frac{d}{2}\right) _{i}\left( d+k-j\right) _{i}}{i!\left( \frac{d}% {2}+k-j+1\right) _{i}\left( \frac{3d}{2}+k+1\right) _{i}}\varepsilon^{2i}. \end{align*} The last sum is a $_{3}F_{2}$ with argument $\varepsilon^{2}.$
Simplify the Gamma terms and note $\left( -k\right) _{j}=\left( -1\right) ^{j}\frac{k!}{\left( k-j\right) !}$ and $\Gamma\left( k-j+1\right) =\left( k-j\right) !.$ Then% \begin{align*} & \frac{\left( -k\right) _{j}}{\left( 1+\frac{d}{2}\right) _{j}}\left( -1\right) ^{j}\frac{\Gamma\left( \frac{d}{2}\right) \Gamma\left( k-j+1\right) \Gamma\left( d+k-j\right) \Gamma\left( \frac{d}% {2}+j+1\right) }{\Gamma\left( \frac{d}{2}+k-j+1\right) \Gamma\left( \frac{3d}{2}+k+1\right) }\\ & =\frac{k!}{\left( 1+\frac{d}{2}\right) _{j}}\frac{\Gamma\left( \frac {d}{2}\right) \Gamma\left( \frac{d}{2}+1\right) \left( \frac{d}% {2}+1\right) _{j}\Gamma\left( d\right) \left( d\right) _{k-j}}% {\Gamma\left( \frac{d}{2}+1\right) \left( \frac{d}{2}+1\right) _{k-j}\Gamma\left( \frac{3d}{2}+1+k\right) }\\ & =\frac{k!\left( d\right) _{k-j}\Gamma\left( \frac{d}{2}\right) \Gamma\left( d\right) }{\left( \frac{d}{2}+1\right) _{k-j}\Gamma\left( \frac{3d}{2}+1+k\right) }. \end{align*} Combine the factors (correctly +/- ??)% \begin{align*} I\left( \varepsilon\right) & =\frac{2\Gamma\left( 1+d+k\right) ^{2}\Gamma\left( \frac{d}{2}\right) \Gamma\left( d\right) k!\varepsilon ^{d}}{d\Gamma\left( \frac{d}{2}\right) ^{3}\Gamma\left( 1+\frac{d}% {2}+k\right) \Gamma\left( 1+k\right) \Gamma\left( \frac{3d}{2}+1+k\right) }\\ & \times\sum_{j=0}^{k}\frac{\left( d\right) _{k-j}}{\left( \frac{d}% {2}+1\right) _{k-j}}\left( 1-\varepsilon^{2}\right) ^{k-j}~_{3}F_{2}\left( % %TCIMACRO{\QATOP{-d/2-j,d/2,d+k-j}{1+d/2+k-j,1+k+3d/2}}% %BeginExpansion \genfrac{}{}{0pt}{}{-d/2-j,d/2,d+k-j}{1+d/2+k-j,1+k+3d/2}% %EndExpansion ;\varepsilon^{2}\right) . \end{align*} The first line simplifies to% [ \frac{\Gamma\left( 1+d+k\right) ^{2}\Gamma\left( 1+d\right) }% {2\Gamma\left( 1+\frac{d}{2}+k\right) \Gamma\left( 1+\frac{d}{2}\right) ^{2}\Gamma\left( 1+\frac{3d}{2}+k\right) }\varepsilon^{d}% ] (apparently agrees with the postulated $k=0$ expression)