# Solve a constrained four-dimensional integration problem for a “magic simplex” probability

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 $$$$\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}$$$$ 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}]]];

• OK, thanks JimB--I noted your comment. I guess my questions tend to be something of a discipline for me--in making a more-or-less permanent, easily retrieved record of results I've obtained and ensuing issues--even if no further progress is made. (I have definitely obtained one collaborator--"student"--through a question mathoverflow.net/questions/325697/… ) – Paul B. Slater Apr 21 '19 at 18:11

3/560 (-4 + 140 ArcTan[4/3] + 35 Sqrt[2] Log[7] -

He has also provided me subsequently with the results of a number of related three-dimensional constrained integration problems for PPT- and "bound-entanglement" probabilities, based mainly on different ranges of the $$q$$ variables.