# Can one test an octonionic interpretation for a conjecture, apparently valid in the complex and quaternionic settings, and proven in the real case?

For the values $\alpha = \frac{1}{2},1, 2$, corresponding to real, complex and quaternionic scenarios, the formulas (https://arxiv.org/abs/1301.6617, eqs. (1)-(3)) $$\label{Hou1} P_1(\alpha) =\Sigma_{i=0}^\infty f(\alpha+i),$$ where $$\label{Hou2} f(\alpha) = P_1(\alpha)-P_1(\alpha +1) = \frac{ q(\alpha) 2^{-4 \alpha -6} \Gamma{(3 \alpha +\frac{5}{2})} \Gamma{(5 \alpha +2})}{3 \Gamma{(\alpha +1)} \Gamma{(2 \alpha +3)} \Gamma{(5 \alpha +\frac{13}{2})}},$$ and $$\label{Hou3} q(\alpha) = 185000 \alpha ^5+779750 \alpha ^4+1289125 \alpha ^3+1042015 \alpha ^2+410694 \alpha +63000 =$$ $$\alpha \bigg(5 \alpha \Big(25 \alpha \big(2 \alpha (740 \alpha +3119)+10313\big)+208403\Big)+410694\bigg)+63000,$$ and (https://arxiv.org/abs/1609.08561, p. 26), with $k=0$, \begin{align*} \label{Hyper1} P_2\left(\alpha\right) & =1-\frac{\alpha\left( 20\alpha+8k+11\right) \Gamma\left( 5\alpha+2k+2\right) \Gamma\left( 3\alpha+k+\frac{3}{2}\right) \Gamma\left( 2\alpha+k+\frac{3}{2}\right) }{2\sqrt{\pi}\Gamma\left( 5\alpha+2k+\frac{7}{2}\right) \Gamma\left( \alpha+k+2\right) \Gamma\left( 4\alpha+k+2\right) }\\ & \times~_{6}F_{5}\left( %TCIMACRO{\QATOP{1,\frac{5}{2}\alpha+k+1,\frac{5}{2}\alpha+k+\frac{3}{2},2\alpha+k+\frac{3}% %{2},3\alpha+k+\frac{3}{2},\frac{5}{2}\alpha+k+\frac{19}{8}}{\alpha+k+2,4\alpha+k+2,\frac{5}% %{2}\alpha+k+\frac{7}{4},\frac{5}{2}\alpha+k+\frac{9}{4},\frac{5}{2}\alpha+k+\frac{11}{8}}}% %BeginExpansion \genfrac{}{}{0pt}{}{1,\frac{5}{2}\alpha+k+1,\frac{5}{2}\alpha+k+\frac{3}{2}% ,2\alpha+k+\frac{3}{2},3\alpha+k+\frac{3}{2},\frac{5}{2}\alpha+k+\frac{19}{8}% }{\alpha+k+2,4\alpha+k+2,\frac{5}{2}\alpha+k+\frac{7}{4},\frac{5}{2}\alpha+k+\frac{9}{4},\frac {5}{2}\alpha+k+\frac{11}{8}}% %EndExpansion ;1\right) \end{align*} both yield to arbitrarily high-precision that $P_1(\frac{1}{2})=P_2(\frac{1}{2})=\frac{29}{64}$, that $P_1(1)=P_2(1)=\frac{8}{33}$, and that $P_1(2)=P_2(2)=\frac{26}{323}$.

Lovas and Andai (https://arxiv.org/abs/1610.01410) have very recently proven that the first of these three cases yields the probability that a random (with respect to Hilbert-Schmidt/Euclidean/flat measure) pair of real quantum bits is separable/unentangled. They also indicated a strategy for proving the second case for pairs of complex (standard) quantum bits. Taking a highly-intensive numerical approach, Fei and Joynt (https://arxiv.org/abs/1409.1993) have found supporting evidence for the three (real, complex and quaternionic) cases.

Uninvestigated, so far, however, is the case $\alpha=4$, for which $P_1(4)=P_2(4)= \frac{44482}{4091349} \approx 0.0108722$ (with, $44482 = 2 \cdot 23 \cdot 967$ and $4091349 = 3 \cdot 29 \cdot 31 \cdot 37 \cdot 41$). This, motivated by random matrix theory, with $\alpha =\frac{\beta}{2}$, with $\beta$ being the usual "Dyson-index", would appear to possibly correspond to some octonionic setting.

So, can one find a framework in which to address the conjecture that $P_1(4)=P_2(4)= \frac{44482}{4091349}$ has an octonionic interpretation?

Let us note that the pairs of real quantum bits constitute the 9-dimensional space of $4 \times 4$ "density" matrices--nonnegative definite, symmetric with real entries and unit trace. The pairs of complex (standard) quantum bits similarly constitute a 15-dimensional space, and the pairs of quaternionic bits, a 27-dimensional space. The pairs of octonionic bits would comprise a 51-dimensional space. (The variable $k$ in the formula for $P_2(\alpha)$, which we set to zero for our discussion, corresponds to a class of ("random induced") measures (https://arxiv.org/abs/quant-ph/0012101), broader than the Hilbert-Schmidt one.)

The two formulas $P_1(\alpha)$ and $P_2(\alpha)$ were developed based solely on analyses of matrices with real and complex (and not quaternionic and octonionic) entries. To be more specific, the ascending moments of determinants (emphasis added) of the 4 x 4 "density" matrices and of their “partial transposes” were computed, and formulas found for them. (These were, then, used in the Mathematica density approximation procedure of Provost [http://www.mathematica-journal.com/issue/v9i4/contents/DensityApproximants/DensityApproximants.pdf], to eventually arrive at the expressions for $P_1(\alpha)$ and $P_2(\alpha)$ .)

The two formulas (Charles Dunkl observed) could be absorbed into one, by regarding the parameter in the complex case to be twice that in the real case (hence the apparent [Dyson-index-like] connection to random matrix theory).

Now, although the calculation of determinants is straightforward with matrices the entries of which are restricted to real and complex values, it becomes more subtle with the quaternions, and, a fortiori, it would seem with the octonions. (Math Reviews has a number of articles, apparently dealing on some level with octonionic determinants.) E. H. Moore (Bull. Amer. Math. Soc. 28 [1922], 161-162) gave a definition in the quaternionic case—and Wikipedia has a brief article, I see, about the “Dieudonne determinant” (“which is a generalization of the determinant of a matrix over division rings and local rings”). Also, the concept of "quasideterminant" (work of Israel Gelfand et al) appears relevant). ("Note also that for octonionic hermitian matrices of size at least 4 no nice notion of determinant is known, while for matrices of size 3 it does exist" S Alesker - Journal of Geometric Analysis, 18 [2008] -[p. 651].)

So, I think the originally stated problem posed above hinges on to what extent the moment formulas Dunkl developed can be “extrapolated” to the octonionic domain. (I note, however, that Fei and Joynt in the cited paper appear to have by-passed the use of determinants, in their quaternionic analysis).

In preparing this “answer”, I found a (rather remarkable) series of June 2012 emails from Dunkl in which he does a highly in-depth [using Maple] analysis of the use of the Moore determinant in the quaternionic case, apparently succeeding in confirming its appropriateness there. Here is part of his treatment (this, of course, deal with the quaternionic scenario, and the octonionic question remains). Dunkl writes more of interest in this series of detailed emails (but I don’t see how to really present his interesting remarks here).

The Maple code Dunkl employed for the application of the Moore determinant to the quaternionic case was:

qm := proc (z1, z2) local zq1, zq2, w1, w2, w3, w4; global qco, iq, jq, kq; zq1 := qco(z1); zq2 := qco(z2); w1 := zq1[1]*zq2[1]-zq1[2]*zq2[2]-zq1[3]*zq2[3 ]-zq1[4]*zq2[4]; w2 := zq1[1]*zq2[2]+zq1[2]*zq2[1]+zq1[3]*zq2[4]-zq1[4]*zq2[3] ; w3 := zq1[1]*zq2[3]+zq1[3]*zq2[1]+zq1[4]*zq2[2]-zq1[2]zq2[4]; w4 := zq1[1] zq2[4]+zq1[4]*zq2[1]+zq1[2]*zq2[3]-zq1[3]*zq2[2]; w1+w2*iq+w3*jq+w4*kq end proc; qconj := proc (f) options operator, arrow; subs({kq = -kq, iq = -iq, jq = -jq} ,f) end proc; qco := proc (f) local fq; global iq, jq, kq, lq; fq := collect(f,[iq, jq, kq]) ; lq[2] := coeff(fq,iq); lq[3] := coeff(fq,jq); lq[4] := coeff(fq,kq); lq[1] := subs(iq = 0,jq = 0,kq = 0,f); [lq[1], lq[2], lq[3], lq[4]] end proc; qdet4x := [[[1, 1], [2, 2], [3, 3], [4, 4], 1], [[4, 4], [3, 3], [1, 2], [2, 1 ], -1], [[4, 4], [2, 3], [3, 2], [1, 1], -1], [[4, 4], [2, 2], [1, 3], [3, 1], -1], [[3, 3], [2, 4], [4, 2], [1, 1], -1], [[3, 4], [4, 3], [2, 2], [1, 1], -1 ], [[3, 4], [4, 3], [1, 2], [2, 1], 1], [[2, 3], [3, 2], [1, 4], [4, 1], 1], [ [2, 4], [4, 2], [1, 3], [3, 1], 1], [[4, 4], [1, 2], [2, 3], [3, 1], -1], [[4, 4], [1, 3], [3, 2], [2, 1], -1], [[3, 3], [1, 4], [4, 2], [2, 1], -1], [[2, 2] , [1, 3], [3, 4], [4, 1], -1], [[2, 2], [1, 4], [4, 3], [3, 1], -1], [[2, 3], [3, 4], [4, 2], [1, 1], -1], [[2, 4], [4, 3], [3, 2], [1, 1], -1], [[1, 2], [2 , 3], [3, 4], [4, 1], 1], [[1, 2], [2, 4], [4, 3], [3, 1], 1], [[1, 3], [3, 2] , [2, 4], [4, 1], 1], [[1, 3], [3, 4], [4, 2], [2, 1], 1], [[1, 4], [4, 2], [2 , 3], [3, 1], 1], [[1, 4], [4, 3], [3, 2], [2, 1], 1], [[3, 3], [2, 2], [1, 4] , [4, 1], -1], [[3, 3], [1, 2], [2, 4], [4, 1], -1]]; qmdet4 := proc (mx) local dt, i, tm, ppq; global qm, qdet4x; dt := 0; for i to 24 do tm := op(i,qdet4x); ppq := qm(mx[tm[1][1],tm[1][2]],qm(mx[tm[2][1],tm[2] [2]],qm(mx[tm[3][1],tm[3][2]],mx[tm[4][1],tm[4][2]]))); dt := dt+tm[5]* simplify(ppq) end do; simplify(dt) end proc;

Also, here is the list of the 24 (4!) factors of the Moore determinant, in order, with the sign:

# this is a list of the 24 factors, in order, with the sign

qdet4x;

[[[1, 1], [2, 2], [3, 3], [4, 4], 1],

    [[4, 4], [3, 3], [1, 2], [2, 1], -1],

[[4, 4], [2, 3], [3, 2], [1, 1], -1],

[[4, 4], [2, 2], [1, 3], [3, 1], -1],

[[3, 3], [2, 4], [4, 2], [1, 1], -1],

[[3, 4], [4, 3], [2, 2], [1, 1], -1],

[[3, 4], [4, 3], [1, 2], [2, 1], 1],

[[2, 3], [3, 2], [1, 4], [4, 1], 1],

[[2, 4], [4, 2], [1, 3], [3, 1], 1],

[[4, 4], [1, 2], [2, 3], [3, 1], -1],

[[4, 4], [1, 3], [3, 2], [2, 1], -1],

[[3, 3], [1, 4], [4, 2], [2, 1], -1],

[[2, 2], [1, 3], [3, 4], [4, 1], -1],

[[2, 2], [1, 4], [4, 3], [3, 1], -1],

[[2, 3], [3, 4], [4, 2], [1, 1], -1],

[[2, 4], [4, 3], [3, 2], [1, 1], -1],

[[1, 2], [2, 3], [3, 4], [4, 1], 1],

[[1, 2], [2, 4], [4, 3], [3, 1], 1],

[[1, 3], [3, 2], [2, 4], [4, 1], 1],

[[1, 3], [3, 4], [4, 2], [2, 1], 1],

[[1, 4], [4, 2], [2, 3], [3, 1], 1],

[[1, 4], [4, 3], [3, 2], [2, 1], 1],

[[3, 3], [2, 2], [1, 4], [4, 1], -1],

[[3, 3], [1, 2], [2, 4], [4, 1], -1]]


# this adds the 24 terms to get the Moore determinant

print(qmdet4);

proc(mx) local dt, i, tm, ppq; global qm, qdet4x; dt := 0; for i to 24 do tm := op(i, qdet4x); ppq := qm(mx[tm[1][1], tm[1][2]], qm( mx[tm[2][1], tm[2][2]], qm( mx[tm[3][1], tm[3][2]], mx[tm[4][1], tm[4][2]]))) ; dt := dt + tm[5]*simplify(ppq) end do; simplify(dt) end proc

# a14,a23,a41,a32

mtqp;

    [    2
[h[1]  , 0 , 0 ,

]
h[2] g2[0] + h[2] g2[1] iq + h[2] g2[2] jq + h[2] g2[3] kq]

[        2
[0 , h[2]  ,

h[1] g1[0] + h[1] g1[1] iq + h[1] g1[2] jq + h[1] g1[3] kq ,

]
0]

[
[0 ,

h[1] g1[0] - h[1] g1[1] iq - h[1] g1[2] jq - h[1] g1[3] kq ,

2        2        2        2       2    ]
g2[0]  + g2[1]  + g2[2]  + g2[3]  + h[3]  , 0]

[
[h[2] g2[0] - h[2] g2[1] iq - h[2] g2[2] jq - h[2] g2[3] kq ,

2        2        2        2       2]
0 , 0 , g1[0]  + g1[1]  + g1[2]  + g1[3]  + h[4] ]


# this is the PT det (Moore formula)

det4;

     2      2       2      2       2      2       2      2


-(-h[2] g2[0] - h[2] g2[1] - h[2] g2[3] - h[2] g2[2]

           2      2       2      2       2      2       2      2
+ h[1]  g1[0]  + h[1]  g1[1]  + h[1]  g1[2]  + h[1]  g1[3]

2     2       2      2       2      2       2      2
+ h[4]  h[1] ) (h[1]  g1[3]  + h[1]  g1[0]  + h[1]  g1[1]

2      2       2     2       2      2       2      2
+ h[1]  g1[2]  - h[2]  h[3]  - h[2]  g2[0]  - h[2]  g2[1]

2      2       2      2
- h[2]  g2[2]  - h[2]  g2[3] )


It appears that we can make meaningful progress in addressing this problem--via numerical computations--but without yet fully resolving a number of issues.

The question as put pertains to $4 \times 4$ (density) matrices--and, in this regard, we seek to extend the quite recent analyses of $2 \times 2$ and $3 \times 3$ "Wishart matrices ($W$) with octonion entries" of Peter Forrester (sec. 3 of https://arxiv.org/pdf/1610.08081.pdf). He employed Cholesky decompositions $W=T^{\dagger} T$. Accordingly, we start with $4 \times 4$ null matrices $T$ and fill their six upper triangular off-diagonal entries with octonions, the eight independent components of each of the six distributed as standard Gaussians.

Next, we fill the four diagonal entries with values that are the square roots of Gamma distribution variates. For the $2 \times 2$ case, Forrester employs $\Gamma[a+1,2]$ and $\Gamma[a+5,2]$. In the $3 \times 3$ instance, he utilizes $\Gamma[a +4 (i-1),2]$, $i=1,2,3$. (The parameter $a$ appears to not need to be fully specified, but must be large enough that the Gamma distributions is well-defined.) For the $4 \times 4$ case, we initially employed $\Gamma[a+3 (i-1),2]$, $i=1,\ldots,4$ (but, admittedly, have no compelling argument for this particular choice).

To further proceed, we relied upon the suite of Mathematica programs made available by Tevian Dray and Corinne A. Manogue in their paper, "Finding octonionic eigenvectors using Mathematica" (Computer Physics Commun 115 [1998], 536-547). This allowed us (using, in particular, their commands "omult[p,q]" and "MMult[X,Y]") to generate random Wishart matrices $W=T^{\dagger} T$.

For each such matrix, we sought to compute its determinant and then test its positivity. For this purpose we attempted to employ "Theorem 5.3. (Laplace expansion)" in the 2010 paper of Jianquan Liao, Jinxun Wang and Xingmin Li in volume 20 of Anal. Theory Appl., entitled "The all-associativity of octonions and its applications"(http://link.springer.com/article/10.1007/s10496-010-0326-2). (We followed the "template" of the Laplace expansion of a $4 \times 4$ matrix by $2 \times 2$ "complementary minors" presented in https://en.wikipedia.org/wiki/Laplace_expansion, and utilized the Mathematica command "Odet[X]" in the Dray-Manogue package for the computation of the $2 \times 2$ minors.)

At this point, we were prepared for our simulation of the Wishart matrices $W=T^{\dagger} T$. At first, we set $a=1$. For one thousand such random matrices, we found (using the Laplace expansion algorithm) all of their determinants to be positive. (Let us note that Forrester in simulating ten thousand $3 \times 3$ Wishart matrices, found 5,500 of them to have negative determinants--but the particular value of $a$ employed was not indicated.)

Now, to address the underlying/central question of the value of the "[Hilbert-Schmidt] separability probability of two-qubit density matrices with octonionic entries", we (again, using the Laplace expansion routine) computed the determinants of the "partial transposes" of the one thousand (positive determinant) Wishart matrices. Of the one thousand partial transposes, 651 had positive determinants, giving us a separability probability of 0.651, orders of magnitude larger than the conjectured value (see initial question) of $P_1(4)=P_2(4) = \frac{44482}{4091349} \approx 0.0108722$.

But now, in our most interesting finding, we ascertained that one could readily "tune" the separability probability by the choice of the parameter $a$. So, for $a=\frac{1}{128}$, we obtained 996 Wishart matrices with positive determinants. Of these, only 13 now had "positive partial transposes (PPT's), giving us a separability probability of $\frac{13}{996} \approx 0.0130522$, much closer to the conjectured value.

So, at this point in time, the main question seems to be does there exist some underlying principle for the choice of the parameter $a$ (as well as what are the appropriate Gamma distributions to employ for setting the diagonal entries of the random upper triangular matrices $T$ in the Cholesky decompositions $W=T^{\dagger} T$ for the random $4 \times 4$ Wishart matrices)? Also, does it suffice to test positivity to only employ determinants, as opposed to the (more computationally problematical) eigenvalues?