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))
begin{equation} label{Hou1}
P_1(alpha) =Sigma_{i=0}^infty f(alpha+i),
end{equation}
where
begin{equation} 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})}},
end{equation}
and
begin{equation} label{Hou3}
q(alpha) = 185000 alpha ^5+779750 alpha ^4+1289125 alpha ^3+1042015 alpha ^2+410694 alpha +63000 =
end{equation}
begin{equation}
alpha bigg(5 alpha Big(25 alpha big(2 alpha (740 alpha
+3119)+10313big)+208403Big)+410694bigg)+63000,
end{equation}
and (https://arxiv.org/abs/1609.08561, p. 26), with $k=0$,
begin{align*} label{Hyper1}
P_2left(alpharight) & =1-frac{alphaleft( 20alpha+8k+11right)
Gammaleft( 5alpha+2k+2right) Gammaleft( 3alpha+k+frac{3}{2}right)
Gammaleft( 2alpha+k+frac{3}{2}right) }{2sqrt{pi}Gammaleft(
5alpha+2k+frac{7}{2}right) Gammaleft( alpha+k+2right) Gammaleft(
4alpha+k+2right) }
& times~_{6}F_{5}left(
%TCIMACRO{QATOP{1,frac{5}{2}alpha+k+1,frac{5}{2}alpha+k+frac{3}{2},2alpha+k+frac{3}%
%{2},3alpha+k+frac{3}{2},frac{5}{2}alpha+k+frac{19}{8}}{alpha+k+2,4alpha+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}%
,2alpha+k+frac{3}{2},3alpha+k+frac{3}{2},frac{5}{2}alpha+k+frac{19}{8}%
}{alpha+k+2,4alpha+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
;1right)
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

here is a PT of a matrix with 4 off-diagonal entries nonzero,

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] ]

the Cholesky factor is [[h1,0,0,g1],[0,h2,g2,0],[0,0,h3,0],[0,0,0,h4]

where g1,g2 are quaternions

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] )