Evaluate 3D- and 5D-constrained integrals for absolute separability probabilities

Question

In a recent posting,

TwoQubits

user JimB, employing a change-of-transformations put forth by N. Tessore, was able to confirm a formula for the "two-qubit absolute separability Hilbert-Schmidt probability" for complex states. This formula had first been expressed in eq. (34) in 2009paper, but with an apparent "typo" there, subsequently corrected in eq. (A2) in

QuasiRandom

Now, in this posting, let us indicate the existence of two companion (clearly more challenging still) problems also discussed in the 2009 paper. One, is again of a 3D nature and the other of a 5D character. For the former, we also have an explicit formula--for which we would desire confirmation and possible simplification--while for the latter, no analytical progress (numerical or symbolic) at all has so far been reported.

---

---

The 3D problem is the quaternionic counterpart to the already successfully addressed complex version. The constraint remains as before (so the transformations of Tessore can, again, be employed to yield an unconstrained problem). However, the powers of two in the integrand are replaced by powers of four, and a new normalization constant is inserted. To now be explicit, again employing $x,y,z$ as the principal variables (rather than subscripted $\lambda$'s for the eigenvalues), the problem takes the form

Integrate[86825246363856000 (x - y)^4 (x - z)^4 (y - z)^4 (-1 + 2 x + y + z)^4 (-1 + x + 2 y + z)^4 (-1 + x + y + 2 z)^4 Boole[x > y && y > z && z > 1 - x - y - z && x - z < 2 Sqrt[y (1 - x - y - z)] && 1 > z > 0 && 1 > y > 0 && 1 > x > 0], {z, 0, 1}, {y, 0, 1}, {x, 0, 1}]

The result of the integration (possibly subject to further simplification) was given in the 2009 paper as

-((13 (s[1] + s[2] + s[3] + s[4] + s[5] + s[6]))/816946343106356485029888)

where, similarly to the two-qubit complex state results, we have the occurrence of $\sqrt{2}$'s and inverse hyperbolic functions,

s[1] = -216449750678398795533760757497856 + 176860737736399592490919645937664 Sqrt[2]


s[2] = 279292548969739228073088142369304501839785 Sqrt[2] Pi


s[3] = -558572941247617043110461841280869072896000 Sqrt[2] ArcCot[Sqrt[2]]




s[4] = 23637916932187025487103667523337320 Sqrt[2] ArcCot[2 Sqrt[2]]


s[5] = -16178155879591789043088455851252390200 Sqrt[2] ArcCot[3 + Sqrt[2]]

and

s[6] = -558589165778586158484606527963549721006600 Sqrt[2] ArcTan[Sqrt[2]].

So, the analytical framework previously successfully employed for the complex states, should be by-and-large applicable, it would seem.

Implementation of the Tessore change-of-variables gives us the equivalent unconstrained integration problem,

Integrate[1/((1 + 2 x)^28 (1 + y)^27) 86825246363856000 (1 + x)^14 (x - y)^4 (1 - 2 z)^4 (y - z)^4 (-1 + y + z)^4 (z + x (-1 - y + z))^4 (-1 + z + x (y + z))^4, {z, 1/2, 1}, {y, z, 2 + 2 Sqrt[1 - z] - z}, {x, y, (4 y + z - 3 y z - z^2 + 2 (1 + y) Sqrt[y - y z])/(-1 + y + z)^2}]

Then, following the structure developed by JimB, one can begin by breaking the problem into five mutually exclusive integrations over x and y , thusly,

integrand = 1/((1 + 2 x)^28 (1 + y)^27) 86825246363856000 (1 + x)^14 (x - 
   y)^4 (1 - 2 z)^4 (y - z)^4 (-1 + y + z)^4 (z + 
   x (-1 - y + z))^4 (-1 + z + x (y + z))^4;


a1 = Integrate[integrand, {y, 1 - z - 2 Sqrt[z - 2 z^2], 1/2 (1 - 2 z)}, {x, 
   1 - y - 2 z, -2 y + z + 2 Sqrt[y - 2 y z]},
  Assumptions -> {1/8 (2 - Sqrt[2]) < z < 1/6}];
a2 = Integrate[integrand, {y, 1/2 (1 - 2 z), (2 - z)/9 + 2/9 Sqrt[1 - z - 2 z^2]}, {x, y, -2 y + z + 2 Sqrt[y - 2 y z]},
  Assumptions -> {1/8 (2 - Sqrt[2]) < z < 1/6}]; 
a3 = Integrate[integrand, {y, z, 1/2 (1 - 2 z)}, {x, 1 - y - 2 z, -2 y + z + 2 Sqrt[y - 2 y z]},
  Assumptions -> {1/6 < z < 1/4}];
a4 = Integrate[integrand, {y, 1/2 (1 - 2 z), (2 - z)/9 + 2/9 Sqrt[1 - z - 2 z^2]}, {x, y, -2 y + z + 2 Sqrt[y - 2 y z]},
  Assumptions -> {1/6 < z < 1/4}]; 
a5 = Integrate[integrand, {y, z, (2 - z)/9 + 2/9 Sqrt[1 - z - 2 z^2]}, {x, 
   y, -2 y + z + 2 Sqrt[y - 2 y z]}];

Then, JimB suggests the use of the rules

rules = {Sqrt[(-1 + 2 z) (-1 + z + 2 Sqrt[(1 - 2 z) z])] -> 1 - 2 z - Sqrt[(1 - 2 z) z],Sqrt[-z (-1 + z + 2 Sqrt[(1 - 2 z) z])] -> Sqrt[(1 - 2 z) z] - z,Sqrt[(1 + z) (2 - z + 2 Sqrt[1 - z - 2 z^2])] -> 1 + z + Sqrt[1 - z - 2 z^2],Sqrt[(-1 + 2 z) (z - 2 (1 + Sqrt[1 - z - 2 z^2]))] -> 1 - 2 z + Sqrt[1 - z - 2 z^2]};

---

---

The 5D problem mentioned at the outset is that of computing the "qubit-qutrit absolute separability Hilbert-Schmidt probability" for complex states. It takes the form

Integrate[8651375866704363561280512000000 (v - w)^2 (v - x)^2 (w - x)^2 (v -  y)^2 (w - y)^2 (x - y)^2 (v - z)^2 (w - z)^2 (x - z)^2 (y - z)^2 (-1 + 2 v + w + x + y + z)^2 (-1 + v + 2 w + x + y + z)^2 (-1 + v + w + 2 x + y + z)^2 (-1 + v + w + x + 2 y + z)^2 (-1 + v + w + x + y + 2 z)^2  Boole[1 > v && v>  w && w > x && x > y && y > z && z > 1 - v - w - x - y - z &&  1 - v - w - x - y - z > 0 && v - z - 2 Sqrt[y (-1 + v + w + x + y + 2 z)] < 0], {z, 0, 1}, {y, 0, 1}, {x, 0, 1}, {w, 0, 1}, {v, 0, 1}]

An extension of the 3D transformation of Tessore to 5D in order to obtain an unconstrained integration would appear appropriate. (Even a numerical evaluation would be of interest.)

In motivating the set of 3D transformations, Tessore had written (possibly helping in the construction of a 5D extension):

"You get the change of variables by first noting that lambda1 + lambda2 + lambda3 < 1, so that lambda1 = x, lambda2 = y (1-x), lambda3 = z (1-y) (1-x) suggests itself. Performing that first change, you see that the limits suggest instead the substitution x -> x/(1+x). Performing that second change, you are then lead to y -> y/(1+x) in the same way. Noting that the resulting range of x at this point is still 0 < x < 1, you can do a substitution x -> x/(1+x) once more to get the change of variables."

Answer 1

This is more a modest insight into the 5D problem than a full solution (the companion 3D problem having been successfully addressed by JimB in the first comment to the twofold question).

Let us break the 5D integration constraint

c = 1 > v && v > w && w > x && x > y && y > z && z > 1 - v - w - x - y - z && 1 - v - w - x - y - z > 0 && v - z - 2 Sqrt[y (-1 + v + w + x + y + 2 z)] < 0

into

c1 = 1 > v && v > w && w > x && x > y && y > z && z > 1 - v - w - x - y - z && 1 - v - w - x - y - z > 0

&&

c2 = v - z - 2 Sqrt[y (-1 + v + w + x + y + 2 z)] < 0 .

Now, the command

h = GenericCylindricalDecomposition[c1, {z, y, x, w, v}][[1]]

yields the answer (Dimensions[h]={2})

(0 < z < 1/6 && ((z < y < 
    1/4 (1 - 
       2 z) && ((y < x < 
        1/3 (1 - y - 2 z) && ((x < w < 1/2 (1 - x - y - 2 z) && 
           1 - w - x - y - 2 z < v < 
            1 - w - x - y - z) || (1/2 (1 - x - y - 2 z) < w < 
            1/2 (1 - x - y - z) && 
           w < v < 1 - w - x - y - z))) || (1/3 (1 - y - 2 z) < 
        x < 1/3 (1 - y - z) && x < w < 1/2 (1 - x - y - z) && 
       w < v < 1 - w - x - y - z))) || (1/4 (1 - 2 z) < y < (
    1 - z)/4 && y < x < 1/3 (1 - y - z) && 
   x < w < 1/2 (1 - x - y - z) && 
   w < v < 1 - w - x - y - z))) || (1/6 < z < 1/5 && z < y < (1 - z)/4 && y < x < 1/3 (1 - y - z) && x < w < 1/2 (1 - x - y - z) && w < v < 1 - w - x - y - z)

The implementation of h[[2]], that is,

1/6 < z < 1/5 && z < y < (1 - z)/4 && y < x < 1/3 (1 - y - z) &&  x < w < 1/2 (1 - x - y - z) && w < v < 1 - w - x - y - z

using an integrand of 1, via the command

Integrate[1, {z, 1/6, 1/5}, {y, z, (1 - z)/4}, {x, y, 1/3 (1 - y - z)}, {w, x, 1/2 (1 - x - y - z)}, {v, w, 1 - w - x - y - z}]

(an unconstrained integration, as Tessore was able to fully achieve in the 3D case) yields

1/111974400  \approx 8.93061*10^-9.

(Using h[[1]] instead yields the considerably larger value of 259/22394880 \approx 0.0000115651, the two results summing to 1/86400, with $86400= 2^7 \cdot 3^3 \cdot 5^2$.)

Now, the modest insight we seek to convey in this "answer" is that the SAME result (1/111974400) is obtained if instead of the constraint c1, we enforce the full constraint c in this integration (so, the vexsome [qubit-qutrit absolute separability] constraint c2 seems irrelevant in this limited context).

In other words, the command

 Integrate[Boole[c2], {z, 1/6, 1/5}, {y, z, (1 - z)/4}, {x, y, 1/3 (1 - y - z)}, {w, x, 1/2 (1 - x - y - z)}, {v, w, 1 - w - x - y - z}]

yields the same result (1/111974400).

So, now instead of an integrand of 1 in the command

    Integrate[1, {z, 1/6, 1/5}, {y, z, (1 - z)/4}, {x, y, 1/3 (1 - y - z)},{w, x, 1/2 (1 - x - y - z)}, {v, w, 1 - w - x - y - z}]

we would like to employ the original integrand

8651375866704363561280512000000 (v - w)^2 (v - x)^2 (w - x)^2 (v - 
y)^2 (w - y)^2 (x - y)^2 (v - z)^2 (w - z)^2 (x - z)^2 (y - 
z)^2 (-1 + 2 v + w + x + y + z)^2 (-1 + v + 2 w + x + y + 
z)^2 (-1 + v + w + 2 x + y + z)^2 (-1 + v + w + x + 2 y + 
z)^2 (-1 + v + w + x + y + 2 z)^2

But such an integration does not seem fully doable--obtaining incomplete beta functions at intermediate steps--with trivariate (x, y, z) results such as

(1/(4 (1 + v1)))(( 2^(-v1 - w1) (2^(2 + v1 + w1) x^(2 + v1 + w1) - (1 - x - y - z)^( 2 + v1 + w1)))/(2 + v1 + w1) + 4 (1 - x - y - z)^( 2 + v1 + w1) (Beta[1/2, 1 + w1, 2 + v1] - Beta[-(x/(-1 + x + y + z)), 1 + w1, 2 + v1])),

where v1 denotes the power of v (w1, x1, y1, z1, similarly) for any of the 152,523 monomials in the expansion of the integrand.

Numerics give us the (infinitesimal-like) probability 1.05552*10^-16--providing a lower bound on the desired absolute separability probability (amounting to some achievement).

I leave this particular integration as a challenge for the community (JimB?)--but I am pessimistic in this regard. Even more so, it would seem, for the additonal required integration--where the c2 constraint is now clearly active--based on the GenericCylindricalDecomposition result h[[1]], rather than h[[2]].