Formulas for Generalized Two-Qubit Separability Probabilities

To begin, we find certain formulas $Q(k,\alpha)= G_1^k(\alpha) G_2^k(\alpha)$, for $k = -1, 0, 1,...,9$. These yield that part of the total separability probability, $P(k,\alpha)$, for generalized (real, complex, quaternionic,\ldots) two-qubit states endowed with random induced measure, for which the determinantal inequality $|\rho^{PT}|>|\rho|$ holds. Here $\rho$ denotes a $4 \times 4$ density matrix, obtained by tracing over the pure states in $4 \times (4 +k)$-dimensions, and $\rho^{PT}$, its partial transpose. Further, $\alpha$ is a Dyson-index-like parameter with $\alpha = 1$ for the standard (15-dimensional) convex set of (complex) two-qubit states. For $k=0$, we obtain the previously reported Hilbert-Schmidt formulas, with (the real case) $Q(0,\frac{1}{2}) = \frac{29}{128}$, (the standard complex case) $Q(0,1)=\frac{4}{33}$, and (the quaternionic case) $Q(0,2)= \frac{13}{323}$---the three simply equalling $ P(0,\alpha)/2$. The factors $G_2^k(\alpha)$ are sums of polynomial-weighted generalized hypergeometric functions $_{p}F_{p-1}$, $p \geq 7$, all with argument $z=\frac{27}{64} =(\frac{3}{4})^3$. We find number-theoretic-based formulas for the upper ($u_{ik}$) and lower ($b_{ik}$) parameter sets of these functions and, then, equivalently express $G_2^k(\alpha)$ in terms of first-order difference equations. Applications of Zeilberger's algorithm yield"concise"forms, parallel to the one obtained previously for $P(0,\alpha) =2 Q(0,\alpha)$. For nonnegative half-integer and integer values of $\alpha$, $Q(k,\alpha)$ has descending roots starting at $k=-\alpha-1$. Then, we (C. Dunkl and I) construct a remarkably compact (hypergeometric) form for $Q(k,\alpha)$ itself. The possibility of an analogous"master"formula for $P(k,\alpha)$ is, then, investigated, and a number of interesting results found.


I. INTRODUCTION
In a previous paper [1], a family of formulas was obtained for the (total) separability probabilities P (k, α) of generalized two-qubit states (N = 4) endowed with Hilbert-Schmidt (k = 0) [2], or more generally, random induced measure [3,4]. In this regard, we note that the natural, rotationally invariant measure on the set of all pure states of a N × K composite system (k = K − N ), induces a unique measure in the space of N × N mixed states [3, eq. (3.6)]. Further, α serves as a Dyson-index-like parameter [5,6], assuming the values 1 2 , 1, 2 for the (N = 4) two-rebit, (standard/complex) two-qubit, and two-quaterbit states, respectively.
The concept itself of a "separability probability", apparently first (implicitly) introduced byŻyczkowski, Horodecki, Sanpera and Lewenstein in their much cited 1998 paper [7], entails computing the ratio of the volume-in terms of a given measure [8]-of the separable quantum states to all quantum states. Here, we first examine a certain component Q(k, α) of P (k, α). This informs us of that portion-equalling simply P (k, α)/2 in the Hilbert-Schmidt (k = 0) case [9]-for which the determinantal inequality |ρ P T | > |ρ| holds, with ρ denoting a 4 × 4 density matrix and ρ P T , its partial transpose. By consequence [10] of the Peres-Horodecki conditions [11,12], a necessary and sufficient condition for separability in this 4 × 4 setting is that |ρ P T | > 0. The nonnegativity condition |ρ| ≥ 0 itself certainly holds, independently of any separability considerations. So, the total separability probability can clearly be expressed as the sum of that part for which |ρ P T | > |ρ| and that for which |ρ| > |ρ P T | ≥ 0. The former quantity will be the one of initial concern here, the ones the formulas Q(k, α) will directly yield.
Here P (k, 1) denotes the total separability probability of the (15-dimensional) standard, complex two-qubit systems endowed with the random induced measure for k = K − 4.

B. Present Analyses
Here, contrastingly ("dually") with respect to the approach indicated in [1], we will find k-specific formulas (k = −1, 0, 1, . . . , 9) as a function of α, that is Q(k, α), for the indicated one (|ρ P T | > |ρ|) of the two component determinantal inequality parts of |ρ P T | > 0. We utilized an exceptionally large number (15,801) number of the first set of moments above in the routine of Provost [17], helping to reveal-to extraordinarily high accuracy-the rational values that the corresponding desired (partial) separability probabilities Q(k, α) strongly appear to assume. Sequences (α = 1, 2, . . . , 30, . . .) of such rational values, then, served as input to the FindSequenceFunction command of Mathematica, which then yielded the initial set of k-specific (hypergeometric-based) formulas for Q(k, α). (This apparently quite powerful [but "black-box"] command of which we have previously and will now make copious use, has been described as attempting "to find a simple function that yields the sequence when given successive integer arguments". It can, it seems, succeed too, at times, for rational-valued inputs.) We, then, decompose Q(k, α) into the product form G k 1 (α)G k 2 (α)
These expressions, in fact, faithfully reproduce the rational-valued (separability probability) sequences that served as input. This fidelity is indicated by numerical calculations to apparently arbitrarily high accuracy (hundreds of digits). (The difference equation results below [sec. V] will provide a basis for our observation as to the rational-valuedness [fractional nature] of these separability probabilities.) In Fig. 1, we show plots of the formulas Q(k, α) obtained over the range α ∈ [1, 10], for Fig. 2, we show a companion plot, exhibiting strongly log-linear-like behavior, for log Q(k, α).
A. Distinguished 7 F 6 function with 2 as an upper parameter in Q(k, α) In each of the eleven k-specific formulas Q(k, α) obtained, there is a distinguished 7 F 6 generalized hypergeometric function, with the ("omnipresent", we will find) argument of (cf. [16]).

The six lower parameters
The lower (bottom) six parameters b ik , i = 1, . . . , 6, of the 7 F 6 function conform for all eleven cases to the simple linear rule, The six entries sum to 6α + 3k + 33 2 .
B. Distinguished 7 F 6 function with 1 as an upper parameter in Q(k, α) Each k-specific formula Q(k, α) we have found also incorporates a second 7 F 6 function (again with argument z = 27 64 , which is, to repeat, invariably the case throughout this paper), having all its thirteen parameters simply equalling 1 less than those in the function just Each of these additional functions possesses, to begin with, the same seven upper parameters (that is, 2, plus those six indicated in (7), (8) and (9)) and the same six lower parameters (6), as in the first 7 F 6 function detailed above (sec. III A). Then, the seven upper parameters are supplemented by from 1 to m k 2's, and the six lower parameters supplemented by from 1 to m k 1's.

11
D. Large α-free terms collapsing to 0 We now point out a rather remarkable property of the formulas for Q(k, α) yielded by the FindSequenceFunction command. If we isolate those (often quite bulky) terms that do not involve any of the m k + 2 hypergeometric functions for each k already described above, we find (to hundreds of digits of accuracy) that they collapse to zero. These terms, typically, do contain hypergeometric functions similar in nature to those described above, but with the crucial difference that the Dyson-index-like parameter α does not occur among their upper and lower parameters. Thus, we are left, after this nullification of terms, with formulas Q(k, α) that are simply sums of m k + 2 polynomial-weighted p F p−1 functions (of α), with p = 7, 7, 8, . . . , 7 + m k .

E. Summary
To reiterate, for each k, our formulas for Q(k, α), all contain a single function of the form There is another distinguished single 7 F 6 function, with all its thirteen parameters being one less. Also there are m k additional functions, i = 1, . . . , m k , with the number of upper 2's running from 2 to m k + 1 and the number of lower 1's, simultaneously running from 1 to m k .
The formulas for Q(k, α) that we have obtained can all be written-we have found-in the product form G k 1 (α)G k 2 (α). The G k 2 (α) factor involves the summation of the hypergeometric functions p F p−1 indicated above, each such function weighted by a polynomial in α, the degrees of the weighting polynomials diminishing as p increases. Let us first discuss the other (hypergeometric-free) factor G k 1 (α), involving ratios of products of Pochhammer symbols.
A. Hypergeometric-function-independent factor G k 1 (α) Some supplementary computations (involving an independent use of the FindSequenceFunction command) indicated that these (hypergeometric-free) factors can be written quite concisely, in terms of the upper and lower parameter sets, setting U ik = where the Pochhammer symbol (rising factorial) is employed. We note that, remarkably, G k 1 (1) = 1-further apparent indication of the special/privileged status of the standard (complex, α = 1) two-qubit states.

Canonical form
In App. A, for k = −1, 0, 1, 2, we show the "canonical form" we have developed for the factors G k 2 (α) (cf. [16, Fig. 3]), the component hypergeometric parts of which we have discussed in sec. III.
It further appears that all the G k 2 (α) factors (k = −1, 0, 1, . . . , 9) (App. B) can be equivalently written as functions that satisfy first-order difference (recurrence) equations of the form where the p's are polynomials in α (cf. [27]). This finding was established by yet another application of the Mathematica FindSequenceFunction command.
That is, we generated-for each value of k under consideration-a sequence (α = 1, 2, . . . , 85) of the rational values yielded by the hypergeometric-based formulas for G k 2 (α), to which the command was then applied. While we have limited ourselves in App. B to displaying our results for k = −1, 0, 1, 2, 3 and 4, we do have the analogous set of results in terms of the hypergeometric functions for the additional instances, k = 5, 6, 7, 8 and 9, and presume that an equivalent set of difference-equation results is constructible (though substantial efforts with k = 5 have not to this point succeeded). The initial points G k 2 (1) in the six difference equations shown are-in the indicated order- Since, as noted above, G k 1 (1) = 1, these are the respective separability probabilities Q(k, 1) themselves. We would like to extend this sequence sufficiently, so that we might be able to establish an underlying rule for it. (However, since the sequence is increasing in value, the Legendre-polynomial densityapproximation procedure of Provost converges more slowly as α increases, so our quest seems somewhat problematical, despite the large number [15,801] to the α-specific values obtained from the so-modified equation to recover the values generated by the original k = −1 difference equation.
A. Polynomial coefficients in difference equations We have for the six (k = −1, 0, 1, 2, 3, 4) cases at hand (App. B) the proportionality where the u ik 's (and b ik 's)-as indicated in sec. III-are themselves functions of α.

The polynomials p k 1 (α)
For all six displayed cases, 14 3. The polynomials p k 0 (α) Further, for all six cases, the polynomial coefficients p k 0 (α)-constituting the inhomogeneous parts of the recurrences-are proportional to the product of a factor of the form and an irreducible polynomial. These irreducible polynomials are, in the indicated order (k = −1, 0, 1), (20) and (for k = 2) The irreducible polynomial for k = 3 is also of degree 7, that is, +698007782α 2 + 471120306α + 134548128.
For k = 4, this auxiliary polynomial p k 0 (α) is now the product of (9+4α) times an irreducible polynomial of degree 7, that is, The coefficients of the highest powers of α in all six irreducible polynomials are factorable into the product of 37 and powers of 2 and 5. 15 In Fig. 3 we show formulas we have generated for the differences between the formulas for Q(k, α) for successive values of k. We note that these are hypergeometric-free. We will find below (51) that these obey the formula Now, as concerns the eleven formulas Q(k, α) (k = −1, 0, 1, . . . , 9) we have obtained for prob(|ρ P T | > |ρ|), which have been the principal focus of the paper, we have computed the ratios of the probability for α = 101 to the probability for α = 100. These ranged from 0.419810 (k = −1) to 0.4204296 (k = 9). Let us note here that z = 27 64 ≈ 0.421875.
We had available α = 1 2 , 1 and 2 computations for k = 1, . . . , 40 for this scenario. We found that, for each of the three values of α, we could construct strongly linear plots-with unit-like slopes between 1.00177 and 1.00297-by taking k times the ratio (R) of the (k + 1) separability probability to the k-th separability probability. (From this, it appears, simply, For values α = k = 1, . . . , 50, we were able to construct a strongly linear plot bysimilarly to the immediate last analysis-taking k = α times the ratio of the (k + 1) = (α + 1) separability probability to the k = α-th separability probability. Now, however, rather than a slope very close to 1, we found a slope near to one-half, that is 0.486882. The (k = α = 0)intercept of the estimated line was 0.894491.
(These results correspond to the variable "dif" in App. C.) Thus, in passing from the (symmetric k = 0) Hilbert-Schmidt setting to the random induced k = 1 scenario, the degree of "conciseness" somewhat diminishes. The polynomials q 0 (α) and q 1 (α) in this pair of formulas are the same as the difference-equation (13) polynomials p 0 0 (α) and p 1 0 (α), given in (19) and (20).

Further, we have
having a root at k = −2.

27
To continue (with a root at k = − 5 2 ), Further, (having a root at k = −3) agreeing with our earlier formulas for all eleven k (as well as k = −2 and 10).
Our formulas give that Q(k, 0) is equal to 1 5 for both k = −2 and 1, and equal to 1 2 for k = 0, . . . , 9. Here, α = 0 presumably corresponds to a classical/nonquantum scenario. Charles Dunkl has observed that for integral values of α, the arguments of the gamma functions in the numerators are of the form 2α + k + 3 2 , 2α + α 2 + k + 3 2 , 3α + k + 3 2 , and in the denominators of the form k + 4n + 1, 2k + 5n + 3 2 . He further noted that the leading (highest power) in the polynomial takes the form α2 5α+2k+1 k α+ α−1 2 −1 . Also, the second leading coefficient (normalizing the leading coefficient of the polynomial to 1) follows the rule Similarly, the so-normalized leading third coefficient takes the form We have been able to generate a considerable number (including k = 1, . . . , 100) of such Q(k, α) formulas, a limited number of which we present in App. D.
Each half-integral α formula contains a gamma function in its numerator with an argument of the form 2 + 5α + 2k and in its denominator a gamma function with an argument of the form 2k + 1 2 (−1) α (−2(−1) α (5α + 2) − i). (In explaining how this formula was obtained, Dunkl stated that the key insights was that Q(k + 1, α) − Q(k, α) factors nicely and that Q(−α − 1, α) = 0.) If we let both α and k be free, and perform the indicated summation in (48), we obtain a hypergeometric-based formula that appears not only to reproduce the formulas in App. D for integer α, but also half-integer and other nonnegative fractional values (such as 1 4 , 2 3 ) of α. Dunkl argued that for k > −α and n = 1, 2, 3, . . .
Taking the limit as n → ∞ The resultant master formula takes the form The value 1 2 from which these terms are subtracted itself has an interesting provenance. It was obtained by conducting the sum indicated in (48), not over j from 0 to α + k as indicated there, but over j from 0 to ∞, that is Q(−α, α) ∞ j=0 H(α, j). (The Q (k, α) formula can then be recovered by subtracting the sum over j from This resulted in the expression (cf. http://math.stackexchange.com/questions/1872364/prove-that-a-certainhypergeometric-function-assumes-either-the-value-frac1) .
For α > 0 this gives us the indicated value of 1 2 . Let us note that for both this 5 F 4 function and the 6 F 5 immediately preceding, the sums of the denominator entries minus the sums of the numerator parameters equal 1 2 -while if these differences had been 1, the two functions could be designated as " 1 2 -balanced" [32]. In the notation of this section (cf. (24)), A. Implications for P (k, α) formula Let us note that for the Hilbert-Schmidt (k = 0) case, apparently [9], 2Q(0, α) = P (0, α), where . Thus, any presumed "master formula" for P (k, α) (sec. XIII), should reduce to 2Q(0, α) for k = 0 (cf. eqs. (25)- (27)). We have been investigating the use of 2Q(k, α) as an initial candidate for P (k, α), then padding out the six upper and five lower entries of the 6 F 5 function with additional pairs of entries, identical for k = 0, but different for k = 0. Then, for k = 0, the initial candidate would be recovered. (The somewhat interesting " 1 2 -balanced" property, mentioned above, or some k-free counterpart of it would, then, be lost.) Initial limited numerical investigations along these lines have been somewhat disappointing, as they appeared to indicate that the best fits would be obtained for pairs of padded entries with equal coefficients of k. Also, fits to values of P (k, α) did not seem to be improved through the padding strategy.
However, another considerably more interesting approach along similarly motivated lines was, then, developed. We mapped the parameter k in the Q(k, α) function to βk, so that for k = 0 the original function would be recovered, no matter the specific value of β. We evaluated the transformed functions by seeing how well they fit the series of (known) eight values P (k, k), k = 5, . . . , 12. For the original β = 1, the figure-of-merit for the fit was 0.7703536. This figure rather dramatically decreases/improves as β increases, reaching a near minimum of 0.0479732 for β = 11 2 (and 0.108008 for β = 5 and 0.153828 for β = 6.) The implications of this phenomenon will be further investigated. Perhaps it might be of value to combine the last two (padding and scaling of k) strategies.

B. Conjectured Identity
In relation to (50), Dunkl formulated the conjecture To avoid zero denominators, it is necessary that α > − 1 8 . For α = 0, the value is 1, while the sum is rational for α = n, n + 1 2 , n = 0, 1, 2 . . .. In response to this conjecture, C. Koutschan wrote: "The 5F4 sum fits into the class of identities that can be done with Zeilberger's algorithm. I attach a Mathematica notebook with some computations. More precisely, using the creative telescoping method, my program finds a linear recurrence equation that is satisfied by the 5F4 sum. It is a trivial calculation to verify that also the right-hand side satisfies the same recurrence. As you remark, both sides give 1 for α = 0. We can conclude that the identity holds for all α in N." However, cases where α is neither an integer or half-integer still require attention. (G. Gasper has commented that the 5 F 4 function is not a special case of the formulas in his paper with M. Rahman [33].)

A. Inspheres
The convex set of two-qubit states possesses an "insphere" of maximum radius. The states within in it are all separable [7,13]. So, one can ask what is the Hilbert-Schmidt separability probability outside of it, presuming the apparent total separability probability of 8  Appendix A in [1] considered the possibility of developing a master formula for the total separability probability P (k, α), that associated with the determinantal inequality (4)). It now clearly seems appropriate to reexamine those results (App. E) in terms of the striking hypergeometric-based formula (sec. XI) we have obtained for the partial separability probability Q(k, α), that associated with the determinantal inequality In the earlier study [1], the formulas took the form of 1 minus terms involving polynomials in k and gamma functions, while above the interesting such terms have been subtracted from 1 2 . So, conjecturally there exists a tightly-related analogue of the results reported in sec. XI for P (k, α). (Dunkl did note the qualitative difference that "the ratio 1 2 −Q(k+1,α) 1 2 −Q(k,a) tends to 1 as k → ∞ but 1−P (k+1,α) 1−P (k,a) tends to 16 27 .") In investigating these matters, we have found that for our set of computed P (k, α), α = 1, . . . , 47, the number and location of the consecutive negative roots (sec X C) are precisely the same (47) as for Q(k, α) (sec. X C). (There strangely appears to be a sole exception to this rule for α = 3, where there are five such roots for Q(k, 3) and six such for P (k, 3), Here, is the equation we have solved to determine-based on [1, App. A]-formulas for P (k, α) for α = 1, . . . , 47. The c's are (nonnegative integer) coefficients we fitted to exact values obtained using the Legendre-polynomial density-approximation routine of Provost [17]. (The first 15,760 of the moments (5) were employed.) .
A. The ratios P (k+1,α)−P (k,α) In App. F, we show a number of formulas we have generated for the differences between the formulas for P (k, α) for successive values of k, in relation to the earlier Q(k, α)-based formulas shown in Fig. 3. (A stark contrast occurs, with the formulas initially yielding rational functions-with equal-degree numerators and denominators-and, then, difference equations.) So, it appears that the quest for a general P (k, α) formula could be successfully addressed by employing the same framework as in the Q(k, α) case, by modifying the H(α, j) function to incorporate the new terms shown in Fig. 3 and their extensions to k, in general.
We see an evident relation between the coefficients of the y[1 + α] terms in the difference equations in App. F and the six hypergeometric upper parameters described in sec. III A 2 in the pattern of two 6's and four 5's. Also, the coefficients of the y[α] terms appear related to the six hypergeometric lower parameters described in sec. III A 1.
Remarkably, when this term was multiplied by the function (which comprises the denominator of the ratio), examples of which are shown in Fig. 3, and formulated in (51), Q(1, α) − Q(0, α) = π2 −4α 3 3α+1 5 −5α−3 (20α + 11)Γ α + 5 6 Γ α + 7 6 Γ(5α + 2) Γ(α)Γ α + 7 10 Γ α + 9 10 Γ α + 11 10 Γ α + 13 10 Γ(2α + 3) , the product simplified to the form 4α(5α+3) 9(α+1) . So, we can consider this term to be the first of two parts of a formula for P (1, α) − P (0, α). Now, in quest of the remaining term, when we formed a new difference equation for just the second part, we obtained a new solution, again naturally breaking into the sum of two parts. Now, the first part-previously given by B. X-states counterpart In App. H we show the analogue of the P (k, α) formulas for the "toy" model of Xstates [28,29]. One feature to be immediately noted is that the arguments of the indicated hypergeometric functions are -1. Another is that for half-integer α's, P (k, α) yields rational values, while P X−states (k, α) yields value of the form 1 minus rational numbers divided by C. Use of consecutive negative roots We have noted that both Q(k, α) and P (k, α) have roots at consecutive negative values of k (sec. X C). If we examine the (limiting) values of P (k, α) for k immediately (one) below the end of the consecutive series, we find that they satisfy the relation (This might serve as a "starting point" analogous to the use ((48), (49)) of Q(−α, α)). For the analogous set of Q(− 1 4 (−1) α ((−1) α (10α + 1) − 1) , α)'s, the real parts appear to be 1 2 for even α and − 1 4 for odd α, with the imaginary parts given by Dunkl has observed that the sequence generated by (60) is really two interspersed sequences, one for odd and one for even values of α. They can be represented as f (2α) = (−1) α and f (2α D. Rules for leading coefficients of the polynomials p α (k) In App. I we show for i = 1, . . . , 10, the first of the rules we have developed for the leading coefficients of the polynomials p α (k) given in the formula above (1) for P (k, α)-having been normalized to monic form (the original leading degree-(4α − 2) coefficient being 2 8α+1 2α−1)! ). (For convenience, we drop this k 4α−2 term, and are left with degree-(4α − 3) polynomials.) We note that these resultant polynomials are of degree 2i. Now, we can make the interesting observation that their leading (highest power) coefficients are given (in descending order) by the rules: Also, C 4 is the product of and Further, C 5 is the product of and −325978342903557i 3 + 2137571940201488i 2 − 8722204904328012i + 13657232612174832.
Continuing, C 6 is the product of and

XIV. CONCLUDING REMARKS
The asymptotic analyses reported here and those in studies of Szarek, Aubrun and Ye [4,36,37] both employ Hilbert-Schmidt and (more generally) random induced measures.
However, contrastingly, we chiefly consider asymptotics as the Dyson-index-like parameter α → ∞ (cf. [38,39]), while they implicitly are concerned with the standard case of α = 1, and large numbers of qubits. Perhaps some relation exists, however, between their high-dimensional findings and the quite limited set of asymptotics we have presented above (secs. VII B, VII C, VIII A 2), pertaining to the dimensional index k → ∞.
A strong, intriguing theme in the analyses presented above has been the repeated occurrence of the interesting constant z = 27 64 = ( 3 4 ) 3 . Let us note that J. Guillera in his article "A new Ramanujan-like series for 1 π 2 ", applying methods related to Zeilberger's algorithm [31],         Exchanges.