Formulas for Generalized Two-Qubit Separability Probabilities

Formulas for Generalized Two-Qubit Separability Probabilities

Paul B. Slater Kavli Institute for Theoretical Physics, University of California, Santa Barbara, CA 93106-4030
July 27, 2019

To begin, we find certain formulas , for . These yield that part of the total separability probability, , for generalized (real, complex, quaternionic,…) two-qubit states endowed with random induced measure, for which the determinantal inequality holds. Here denotes a density matrix, obtained by tracing over the pure states in -dimensions, and , its partial transpose. Further, is a Dyson-index-like parameter with for the standard (15-dimensional) convex set of (complex) two-qubit states. For , we obtain the previously reported Hilbert-Schmidt formulas, with (the real case) , (the standard complex case) , and (the quaternionic case) —the three simply equalling . The factors are sums of polynomial-weighted generalized hypergeometric functions , , all with argument . We find number-theoretic-based formulas for the upper () and lower () parameter sets of these functions and, then, equivalently express in terms of first-order difference equations. Applications of Zeilberger’s algorithm yield “concise” forms of and , parallel to the one obtained previously (J. Phys. A, 46 [2013], 445302) for . For nonnegative half-integer and integer values of , (as well as ) has descending roots starting at . Then, we (C. Dunkl and I) construct a remarkably compact (hypergeometric) form for itself. The possibility of an analogous “master” formula for is, then, investigated, and a number of interesting results found.

quantum systems, entanglement probability distribution moments, probability distribution approximation, Peres-Horodecki conditions, partial transpose, determinant of partial transpose, two qubits, two rebits, induced measures, Hilbert-Schmidt measure, moments, separability probabilities, determinantal moments, inverse problems, random matrix theory, generalized two-qubit systems, hypergeometric functions, difference equations, large beta asymptotics, Zeilberger’s algorithm, creative telescoping, Dyson indices, FindSequenceFunction command
Valid PACS 03.67.Mn, 02.30.Zz, 02.50.Cw, 02.40.Ft, 03.65.-w

I Introduction

In a previous paper Slater and Dunkl (2015a), a family of formulas was obtained for the (total) separability probabilities of generalized two-qubit states () endowed with Hilbert-Schmidt () Życzkowski and Sommers (2003), or more generally, random induced measure Życzkowski and Sommers (2001); Aubrun et al. (2014). In this regard, we note that the natural, rotationally invariant measure on the set of all pure states of a composite system (), induces a unique measure in the space of mixed states (Życzkowski and Sommers, 2001, eq. (3.6)). Further, serves as a Dyson-index-like parameter Dyson (1970); Dumitriu et al. (2007), assuming the values for the () 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 Życzkowski et al. (1998), entails computing the ratio of the volume–in terms of a given measure Petz and Sudár (1996)–of the separable quantum states to all quantum states. Here, we first examine a certain component of . This informs us of that portion–equalling simply in the Hilbert-Schmidt () case Slater and Dunkl (2015b)–for which the determinantal inequality holds, with denoting a density matrix and , its partial transpose. By consequence Augusiak et al. (2008) of the Peres-Horodecki conditions Peres (1996); Horodecki et al. (1996), a necessary and sufficient condition for separability in this setting is that . The nonnegativity condition 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 and that for which . The former quantity will be the one of initial concern here, the ones the formulas will directly yield.

The complementary quantity, that for which can, in the most basic cases of interest, be readily obtained from the total separability probability formulas reported in Slater and Dunkl (2015a), which took the form


where for integral and half-integral ,


Here, for integral , is a polynomial of degree with leading coefficient

In Slater and Dunkl (2015a), certain -specific formulas ( and ) had been derived (and we have since continued the integral series to ). Most notably (Slater and Dunkl, 2015a, eq. (3)),


Here denotes the total separability probability of the (15-dimensional) standard, complex two-qubit systems endowed with the random induced measure for . Further, in the two-quater[nionic]bit setting (Slater and Dunkl, 2015a, eq. (4)),


Also, for the two-re[al]bit scenario (Slater and Dunkl, 2015a, eq. (5)),


Tables 1, 2 and 3 in Slater and Dunkl (2015a) reported for , the, in general, rather simple fractional separability probabilities yielded by these three formulas.

By way of example, we first note that formula (2) yields . Then, since we will find from our analyses below, that , we can readily deduce that the corresponding (complementary) separability probability corresponding to the inequalities , for this scenario is equal to .

Let us further observe that for the Hilbert-Schmidt () case, strong evidence has been presented Slater and Dunkl (2015b) that for the two-rebit, two-qubit and two-quaterbit cases, the apparent total separability probabilities of and , respectively, are equally divided between the two forms of determinantal inequalities (cf. Szarek et al. (2006). Lovas and Andai have recently formally proven this two-rebit result and presented an integral formula they hope to similarly yield the two-qubit proportion Lovas and Andai (2016). (These “half-probabilities”, remarkably, are also the corresponding separability probabilities of the minimally degenerate states Szarek et al. (2006), those for which has a zero eigenvalue.) For , however, our analyses will indicate that equal splitting is not, in fact, the case. Greater separability probability is associated with the inequality than . Thus, in the instance just discussed, we do have . (On the other hand, if , then necessarily , so all the total separability probability must, it is clear, be assigned to the component. That is, .) Observations of this nature should help in the further understanding of the intricate geometry of the generalized two-qubit states endowed with random induced measure (cf. Gamel (2016)).

Ii Procedures

ii.1 Previous Analyses

To obtain the new formulas to be presented here for the separability probability amounts for which holds, we first employed–as in our prior studies Slater and Dunkl (2012); Slater (2013); Slater and Dunkl (2015b); Slater and Dunkl (2015a)–the Legendre-polynomial-based probability density approximation (Mathematica-implemented) algorithm of Provost Provost (2005) (cf. Askey et al. (1982)). In this regard, we utilized the previously-obtained determinantal moment formula (Slater and Dunkl, 2015a, eq. (6)) (Slater and Dunkl, 2015b, sec. II) (cf. Bartkiewicz et al. (2015))

(where the variable has the same sense as indicated above, in equalling , and the bracket notation indicates averaging with respect to the random induced measure). Here, , where the Pochhammer (rising factorial) notation is employed.

On the other hand in Slater and Dunkl (2015a), a second companion moment formula (Slater and Dunkl, 2012, sec. X.D.6)


had been utilized for density-approximation purposes with the routine of Provost, with the objective of finding the total separability probabilities , associated with the Peres-Horodecki-based inequality . (These moment formulas had been developed in Slater and Dunkl (2012), based on calculations solely for the two-rebit [] and two-qubit [] cases. However, they do appear, as well, remarkably, to apply to the two-quater[nionic]bit [] case, as reported by Fei and Joynt in a highly computationally intensive Monte Carlo study Fei and Joynt (). No explicit formal extension of the Peres-Horodecki positive-partial-transposition separability conditions Peres (1996); Horodecki et al. (1996) to two-quaterbit systems seems to have been developed, however [cf. Caves et al. (2001); Aslaksen (1996); Peres (1979)]. The value corresponds, presumably it would seem, to an octonionic setting Najarbashi et al. (2016); Forrester (2016).)

ii.2 Present Analyses

Here, contrastingly (“dually”) with respect to the approach indicated in Slater and Dunkl (2015a), we will find -specific formulas () as a function of , that is , for the indicated one () of the two component determinantal inequality parts of . We utilized an exceptionally large number (15,801) number of the first set of moments above in the routine of Provost Provost (2005), helping to reveal–to extraordinarily high accuracy–the rational values that the corresponding desired (partial) separability probabilities strongly appear to assume. Sequences () of such rational values, then, served as input to the FindSequenceFunction command of Mathematica, which then yielded the initial set of -specific (hypergeometric-based) formulas for . (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, and perhaps even ones of a symbolic nature.) We, then, decompose into the product form

Iii Common features of the -specific formulas

For each , the FindSequenceFunction command yielded what we can consider as a large, rather cumbersome (several-page) formula, which we denote by . 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 obtained over the range , for . For fixed , we have , if . In Fig. 2, we show a companion plot, exhibiting strongly log-linear-like behavior, for .

Figure 1: Plots of the separability probability formulas over the range , for . For fixed , we have , if .
Figure 2: Plots of over the range , for . For fixed , , if .

iii.1 Distinguished function with 2 as an upper parameter in

In each of the eleven -specific formulas obtained, there is a distinguished generalized hypergeometric function, with the (“omnipresent”, we will find) argument of (cf. Guillera (2011) (Koepf, 2014, Ex. 8.6, p. 159)), having 2 as one of the seven upper parameters (cf. Slater (2013)).

iii.1.1 The six lower parameters

The lower (bottom) six parameters , , of the function conform for all eleven cases to the simple linear rule,


The six entries sum to .

iii.1.2 The six upper parameters

The six upper parameters (aside from the seventh -invariant constant of 2 already indicated), , can be broken into one set of two (the numerical parts summing to integers), incorporating consecutive fractions having 6’s in their denominators, and one set of four (the numerical parts also summing to integers), incorporating consecutive fractions having 5’s in their denominators. For the set of two, the smaller of the two upper entries abides by the rule


where the (integer-valued) floor function is employed, while the larger entry is given by


For integral values of , the same values of and are yielded by the interpolating functions,



For , for illustrative purposes, application of the two rules yields , and for , we have . (We have noted that is an integer. The sequence of these integers–for arbitrary integer or half-integer values of –is found in the On-Line Encyclopedia of Integer Sequences [] as A004523 [“Two even followed by one odd”] and as A232007 [“Maximal number of moves needed to reach every square by a knight from a fixed position on an n X n chessboard, or -1 if it is not possible to reach every square”].)

For the complementary set of four upper parameters of the function, the entries in order of increasing magnitude are expressible as



For , for illustrative purposes, application of these four rules yields , and for , we have . For arbitrary , the sum of the four terms under discussion minus is an integer, namely, . Further, let us note that for integral values of , has values


iii.2 Distinguished function with 1 as an upper parameter in

Each -specific formula we have found also incorporates a second function (again with argument , 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 described. (A basic transformation exists [consulting the HYP manual of C. Krattenthaler, available at], allowing one to convert the thirteen [twelve -dependent parameters, plus 1] of this function [that is, add 1 to each of them] to those thirteen of the first distinguished previously described, plus other terms.)

iii.3 The remaining functions, , in .

Now, in addition to the two distinguished functions just presented, there are more hypergeometric functions , , for each , where


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 function detailed above (sec. III.1). Then, the seven upper parameters are supplemented by from 1 to 2’s, and the six lower parameters supplemented by from 1 to 1’s.

iii.4 Large -free terms collapsing to 0

We now point out a rather remarkable property of the formulas for yielded by the FindSequenceFunction command. If we isolate those (often quite bulky) terms that do not involve any of the hypergeometric functions for each 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 that are simply sums of polynomial-weighted functions (of ), with .

iii.5 Summary

To reiterate, for each , our formulas for , all contain a single function of the form


There is another distinguished single function, with all its thirteen parameters being one less. Also there are additional functions, ,

with the number of upper 2’s running from 2 to and the number of lower 1’s, simultaneously running from 1 to .

Iv Decomposition of into the product

The formulas for that we have obtained can all be written–we have found–in the product form . The factor involves the summation of the hypergeometric functions indicated above, each such function weighted by a polynomial in , the degrees of the weighting polynomials diminishing as increases. Let us first discuss the other (hypergeometric-free) factor , involving ratios of products of Pochhammer symbols.

iv.1 Hypergeometric-function-independent factor

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 , as


where the Pochhammer symbol (rising factorial) is employed. We note that, remarkably, –further apparent indication of the special/privileged status of the standard (complex, ) two-qubit states.

iv.2 Hypergeometric-function-dependent factor

iv.2.1 Canonical form

In App. A, for , we show the “canonical form” we have developed for the factors (cf. (Slater, 2013, Fig. 3)), the component hypergeometric parts of which we have discussed in sec. III.

V Difference equation formulas for

It further appears that all the factors () (App. B) can be equivalently written as functions that satisfy first-order difference (recurrence) equations of the form


where the ’s are polynomials in (cf. Petkovšek (1992)). This finding was established by yet another application of the Mathematica FindSequenceFunction command.

That is, we generated–for each value of under consideration–a sequence () of the rational values yielded by the hypergeometric-based formulas for , to which the command was then applied. While we have limited ourselves in App. B to displaying our results for and 4, we do have the analogous set of results in terms of the hypergeometric functions for the additional instances, and 9, and presume that an equivalent set of difference-equation results is constructible (though substantial efforts with have not to this point succeeded). The initial points in the six difference equations shown are–in the indicated order–. The next five members of this monotonically-increasing sequence are . Since, as noted above, , these are the respective separability probabilities 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 density-approximation procedure of Provost converges more slowly as increases, so our quest seems somewhat problematical, despite the large number [15,801] of moments incorporated [cf. (Slater and Dunkl, 2015a, App. II)].)

If in the difference equation for , we replace the term by , then we can add


to the -specific values obtained from the so-modified equation to recover the values generated by the original difference equation.

v.1 Polynomial coefficients in difference equations

v.1.1 The polynomials

We have for the six () cases at hand (App. B) the proportionality relation


where the ’s (and ’s)–as indicated in sec. III–are themselves functions of .

v.1.2 The polynomials

For all six displayed cases,


v.1.3 The polynomials

Further, for all six cases, the polynomial coefficients –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 (),


and (for )


The irreducible polynomial for is also of degree 7, that is,


For , this auxiliary polynomial is now the product of 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.

Vi Hypergeometric-Free Formulas for

In App. C we show formulas we have generated for the differences between the formulas for for successive values of . We note that these are hypergeometric-free. We will find below (51) that these obey the formula


Vii Partial separability probability asymptotics

vii.1 -specific formulas

Now, as concerns the eleven formulas () we have obtained for , which have been the principal focus of the paper, we have computed the ratios of the probability for to the probability for . These ranged from 0.419810 () to 0.4204296 (). Let us note here that .

vii.2 -specific formulas

We had available and 2 computations for 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 times the ratio () of the separability probability to the -th separability probability. (From this, it appears, simply, that , as .)

vii.3 “Diagonal” formulas

For values , we were able to construct a strongly linear plot by–similarly to the immediate last analysis–taking times the ratio of the separability probability to the -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 ()-intercept of the estimated line was 0.894491.

Viii Total separability probability formulas

Efforts of our to conduct parallel sets of (-specific) analyses to those reported above for the total separability probabilities , corresponding to , rather than for that component part of the probabilities satisfying the determinantal inequality had been unsuccessful, in the following sense. We had computed what appeared to be appropriate sequences of rational values for and for , but the Mathematica FindSequenceFunction did not yield any underlying governing rules. (This can be contrasted with the results in Slater and Dunkl (2015a), where such successes were reported in obtaining -specific [] formulas [ and ], including (2)-(4) above. However, we do eventually succeed in characterizing the nature of these two () sequences [cf. sec. G].)

In Fig. 3, we plot the logs of these seventy-four total separability probabilities (based on ). A least-squares linear fit to these points is , while in Fig. 4, we show (based on ) the counterpart, with an analogous fit of . (We note that .) Although the slopes of these two linear fits are quite close, the -intercepts themselves are of different sign. The predicted probabilities at , the first of the fitted points, are 0.289019 and 0.602955, respectively. In statistical parlance, the “coefficients of determination” or for the two linear fits to the log-plots are both greater than 0.99995. Further, sampling at , we obtained an estimated, again, very-well fitting line of .

Figure 3: Plot of logs of total separability probability () for random induced measure with . A least-squares linear fit to these 74 points is .
Figure 4: Plot of logs of total separability probability () for random induced measure with . A least-squares linear fit to these 124 points is . We note that .

viii.1 Total separability probability asymptotics

viii.1.1 -specific formulas

C. Dunkl, on the basis of our analysis just above (and its companions), did advance the bold and (certainly, in our overall analytical context) elegant hypothesis of a -invariant () slope equal to , which does seem quite consistent with the numerical properties we have observed (that is, with the direction in which the estimates of the slope tend as the number of points sampled increase).

As further support, we obtained for a analysis, a slope estimate of -0.864025, again converging in the direction of . (Let us remark, regarding the generalized two-qubit version of the [simpler, lower-dimensional] X-states model Mendonça et al. (2014); Dunkl and Slater (2015); Bartkiewicz et al. (2015), that it has been shown that the slope of a [now, log-log] plot of vs. tends to , as .)

viii.1.2 -specific formulas

These interesting observations led us to reexamine, for their asymptotic properties, the “dual” formulas (2)-(4), given above, and previously reported in Slater and Dunkl (2015a). We now find–through analytic means–that for each of and , that as , the ratio of the logarithm of the -st separability probability to the logarithm of the -th separability probability is (cf. (Chu and Zhang, 2014, sec. 7)). (Presumably, the pattern continues for larger , but the required computations have, so far, proved too challenging.)

For example, for , we have for the two-rebit total separability probability , as a function of , the formula (4) given above. In Fig. 5, we show a plot of vs. . The slope of a least-squares-fitted line based on the 200 points is -0.523280, while . (As we increase from , but hold the number of points constant at 200, the approximation of the slope to this value slowly weakens.)

Figure 5: Plot of vs. . The slope of a least-squares-fitted line is -0.523280, while .

Ix “Concise formulas” for

Let us remind the reader of the interesting “concise” (Hilbert-Schmidt []) generalized two-qubit result–applying Zeilberger’s (“creative telescoping”) algorithm Zeilberger (1990)–of Qing-Hu Hou, reported in (Slater, 2013, eqs. (1)-(3)). This–in our present notation–takes the form






We divide the originally reported formula by one-half Slater and Dunkl (2015b), since we have moved here from the () Hilbert-Schmidt original scenario to its counterpart. Using our earlier results above, Hou has further been able to construct the analogue of the “concise formula” above (a Maple worksheet of his is presented in App. D [ cf. (Slater, 2013, Figs. 5, 6)). That is,






(These results correspond to the variable “dif” in App. D.) Thus, in passing from the (symmetric ) Hilbert-Schmidt setting to the random induced scenario, the degree of “conciseness” somewhat diminishes. The polynomials and in this pair of formulas are the same as the difference-equation (13) polynomials and , given in (19) and (20).

At this point in our research, we were able to employ the Mathematica-based HolonomicFunctions package of Christoph Koutschan of the Research Institute for Symbolic Computation (RISC) of Johannes Kepler University. With it, we were readily able to derive the result




We see that the polynomial above is, in expanded form, the same as given in (18).

For the standard trio of Dyson-indices and 2, this formula for yields