Equation (7) in the 2012 paper, "Complementarity Reveals Bound Entanglement of Two Twisted Photons" of B. C. Hiesmayr and W. Löffler for a state $\rho_d$ in the "magic simplex" of Bell states \begin{equation} \rho_d= \frac{q_4 (1-\delta (d-3)) \sum _{z=2}^{d-2} \left(\sum _{i=0}^{d-1} P_{i,z}\right)}{d}+\frac{q_2 \sum _{i=1}^{d-1} P_{i,0}}{(d-1) (d+1)}+\frac{q_3 \sum _{i=0}^{d-1} P_{i,1}}{d}+\frac{\left(-\frac{q_1}{d^2-d-1}-\frac{q_2}{d+1}-(d-3) q_4-q_3+1\right) \text{IdentityMatrix}\left[d^2\right]}{d^2}+\frac{q_1 P_{0,0}}{d^2-d-1} \end{equation} yields for certain values of the $q_i$'s, "for $d=3$ the one-parameter Horodecki-state, the first found bound entangled state".
More generally, for the case $d=3$, the constraint requiring that the partial transpose (obtained by transposing in place the nine $3 \times 3$ blocks) of the density matrix $\rho_3$ be positive definite takes the form
constraint3=q1>0&&q2>0&&q3>0&&4 q1+5 q2+20 q3<20&&512 q1^2+80 q1 (8-11 q2+4 q3)+25 (5 q2^2+16 q2 (2+q3)+64 (-1+q3) (1+2 q3))<0
The command
Integrate[Boole[constraint3],{q1,0,5},{q2,0,4},{q3,0,1}]/(10/3)
then, interestingly, yields the Hilbert-Schmidt "PPT-probability" that the partial transpose of $\rho_3$ is positive definite,
(1/13720)(-4312 + 5145 \[Pi] + 2240 Sqrt[7] ArcCos[11/(8 Sqrt[2])] - 5160 Sqrt[7] ArcSin[(5 Sqrt[7])/16] - 6860 ArcTan[7] + 6280 Sqrt[7] ArcTan[(5 Sqrt[7])/9])
which is approximately 0.461554. (This result was posted as a comment to my earlier query https://quantumcomputing.stackexchange.com/posts/5943/edit .)
Now, I would like to similarly solve the still more formidable $d=4$ problem. Then, the constraint (found by enforcement of the positive definiteness of the sixteen nested minors of both $\rho_4$ and its partial transpose $\rho_4^{PT}$) takes the form
constraint4=q1>0&&q2>0&&q3>0&&q4>0&&5 q1+11 (q2+5 (q3+q4))<55&&3375 q1^2+121 (7 q2^2+90 q2 (1+q3-q4)+225 (1+3 q3-q4) (-1+q3+q4))<330 q1 (19 q2-15 (1+q3-q4))&&(45 q1+11 (15-7 q2-15 q3+45 q4)) (75 q1-11 (15+q2-15 q3+45 q4))<0
Enforcement of the command
Integrate[Boole[constraint4],{q3,0,1},{q2,0,5},{q1,0,11},{q4,0,1}]/(55/24)
would, then, yield the corresponding Hilbert-Schmidt PPT-probability. (The particular ordering of the four variables was suggested by the GenericCylindricalDecomposition command for the 24 possible orderings, but, of course, variations can be investigated.)
Presently, using simply the free form of the WolframCloud, my various attempts to perform the integration--by one approach or another--get timed out. In any case, the problem may be too formidable, by any means. (Perhaps some transformations of variables could be effective.)
Given such PPT-probabilities, the next question that would arise--of a nature that has never really been significantly addressed--is how the probabilities are divided between "bound entangled" and "separable" states (see Fig. 3 of the cited Hiesmayr/Löffler paper).
This code can be employed to generate $\rho_4$
d = 4; W[k_, l_] := Sum[Exp[2 Pi I k n/d] Outer[Times, S[n], S[Mod[n + l, 4]]], {n, 0, d - 1}]; S[0] = {1, 0, 0, 0}; S[1] = {0, 1, 0, 0}; S[2] = {0, 0, 1, 0}; S[3] = {0, 0, 0, 1}; Omega[0, 0] = (1/Sqrt[4]) Sum[ TensorProduct[S[s], S[s]], {s, 0, d - 1}]; Do[ Omega1[k, l] = ArrayReshape[ TensorProduct[W[k, l], IdentityMatrix[4] Omega[0, 0]], {16, 16}]/ Sqrt[4], {k, 0, d - 1}, {l, 0, d - 1}]; Do[ P[k, l] = Outer[Times, Omega1[k, l].ConjugateTranspose[Omega1[k, l]]], {k, 0, d - 1}, {l, 0, d - 1}]; den = Sum[c[k, l] P[k, l], {k, 0, d - 1}, {l, 0, d - 1}];rho[d_] := (1 - q1/(d^2 - (d + 1)) - q2/(d + 1) - q3 - (d - 3) q4) IdentityMatrix[d^2]/d^2 + q1 P[0, 0]/(d^2 - (d + 1)) + q2/((d + 1) (d - 1)) Sum[P[i, 0], {i, 1, d - 1}] + (q3/d) Sum[ P[i, 1], {i, 0, d - 1}]+(q4/d) Sum[Sum[P[i,z],{i,0,d-1}],{z,2,d-2}];rho[4]
The partial transpose of $\rho_4$ is obtained by
ArrayFlatten[Transpose[Partition[rho[4], {4,4}]]];
