# Find the probability (relative volume) of a certain 4-ball with respect to Hilbert-Schmidt measure

Let us consider the set of points {x,y,z,1-x-y-z} and impose the strict ordering constraint

1 > x && x > y && y > z && z > (1 - x - y - z) && (1 - x - y - z) > 0 .

We are interested in the formula for the relative volume/probability ($$\approx 0.0483353$$) of the 4-ball,

x^2 + y^2 + z^2 + (1 - x - y - z)^2 < 3/8,

subject to this ordering constraint, with respect to the Hilbert-Schmidt measure HSmeasure,

9081072000 (x - y)^2 (x - z)^2 (y - z)^2 (-1 + 2 x + y + z)^2 (-1 + x + 2 y + z)^2 (-1 + x + y + 2 z)^2  .

(That is, we want to integrate this measure over the set defined by the intersection of these two constraints.)

If we replace 3/8 by 1/3 in the problem, we have computed the associated Hilbert-Schmidt probability to be $$$$\frac{35 \pi }{23328 \sqrt{3}} \approx 0.00272132 .$$$$ (Let us note that along with $$35 =5 \cdot 7$$, we have $$23328 =2^5 \cdot 3^6$$.)

The larger ball circumscribes and the smaller ball inscribes the set of "two-qubit absolutely separable states" Adhikari, for which the relevant defining constraint, rather than involving a sum of squares, is

x - z < 2 Sqrt[y (1 - x - y - z)]  .

The Hilbert-Schmidt probability of this set has been shown by user JimB in his answer AbsSep to be

29902415923/497664 - 50274109/(512 Sqrt[2]) - (3072529845 π)/(32768 Sqrt[2]) +(1024176615 ArcCos[1/3])/(4096 Sqrt[2])

$$\approx 0.00365826$$.

We suspect/conjecture that, as for the smaller inscribed ball, the sought answer for the larger ball is also a multiple of $$\pi$$ (with possibly also interestingly factorizable denominator and numerator coefficients).

To define probabilities above, we compute volumes with respect to the volume of all two-qubit states (absolutely separable and otherwise). That is, the integral--enforcing only the strict ordering constraint--of the indicated Hilbert-Schmidt measure over the entirety of two-qubit states is 1.

Perhaps the 3D constrained integration can be converted into an unconstrained problem, as Tessore was able to accomplish in AbsSep.

I have not yet investigated if a transformation to hyperspherical coordinates might prove useful. `

Of course, one could pose a still more challenging problem in which the radii squared (that is, 1/3 and 3/8) is simply replaced by a variable.

It appears that the "3/8-problem" posed is more challenging than the "1/3-problem" previously solved--in ways I intend to detail.

(35 Sqrt[1/3 (2692167889921345 - 919847607929856 Sqrt[6])] \[Pi])/27518828544

\approx 0.0483353 .

As speculated in the question--based on the previously derived result of

(35 \[Pi])/(23328 Sqrt[3])

\approx 0.00272132

x^2 + y^2 + (1 - x - y - z)^2 + z^2 < 1/3

was used, rather than

x^2 + y^2 + (1 - x - y - z)^2 + z^2 < 3/8

the new answer does have a factor of $$\pi$$ and an interestingly-factorizable denominator (27518828544) equal to $$2^{22} \cdot 3^8$$.

To derive the result, I employed the command

m = GenericCylindricalDecomposition[1 > x && x > y && y > z && z > (1 - x - y - z) && (1 - x - y - z) > 0&&x^2 + y^2 + z^2 + (1 - x - y - z)^2 < 3/8, {z,x,y}][[1]];

yielding

(1/8 (2 - Sqrt[2]) < z <
1/24 (6 - Sqrt[
6]) && ((1/2 (1 - 2 z) - 1/4 Sqrt[-1 + 16 z - 32 z^2] < y <
1/2 (1 - 2 z) &&
1 - y - 2 z < x <
1/2 (1 - y - z) +
1/4 Sqrt[-1 + 8 y - 12 y^2 + 8 z - 8 y z - 12 z^2]) || (1/
2 (1 - 2 z) < y < (1 - z)/3 + 1/12 Sqrt[1 + 16 z - 32 z^2] &&
y < x <
1/2 (1 - y - z) +
1/4 Sqrt[-1 + 8 y - 12 y^2 + 8 z - 8 y z - 12 z^2]))) || (1/
24 (6 - Sqrt[6]) < z < 1/
6 && ((z < y < 1/2 (1 - 2 z) &&
1 - y - 2 z < x <
1/2 (1 - y - z) +
1/4 Sqrt[-1 + 8 y - 12 y^2 + 8 z - 8 y z - 12 z^2]) || (1/
2 (1 - 2 z) < y < (1 - z)/3 + 1/12 Sqrt[1 + 16 z - 32 z^2] &&
y < x <
1/2 (1 - y - z) +
1/4 Sqrt[-1 + 8 y - 12 y^2 + 8 z - 8 y z - 12 z^2]))) || (1/
6 < z < 1/
4 && ((z < y < 1/2 (1 - 2 z) &&
1 - y - 2 z < x <
1/2 (1 - y - z) +
1/4 Sqrt[-1 + 8 y - 12 y^2 + 8 z - 8 y z - 12 z^2]) || (1/
2 (1 - 2 z) < y < (1 - z)/2 - 1/4 Sqrt[-1 + 8 z - 12 z^2] &&
y < x < 1/2 (1 - y - z) +
1/4 Sqrt[-1 + 8 y - 12 y^2 + 8 z - 8 y z - 12 z^2]) || ((
1 - z)/2 - 1/4 Sqrt[-1 + 8 z - 12 z^2] < y < (1 - z)/2 &&
y < x < 1 - y - z))) || (1/4 < z < 1/3 && z < y < (1 - z)/2 && y < x < 1 - y - z)

I performed the first two integrations (over x, followed by y).

This left me with 270 terms in z, with four ranges of z over which to be integrated.

Many of these terms required use of the "denestSqrt" function, presented by Carl Woll in his answer to denestSqrt .

Among them were integrations of the form

Integrate[(6218205593600 z^i ArcTan[(2 - 2 z - 3 Sqrt[-1 + 8 z - 12 z^2])/(Sqrt[3] Sqrt[3 - 8 z^2 + 4 Sqrt[-1 + 8 z - 12 z^2] - 4 z Sqrt[-1 + 8 z - 12 z^2]])])/(729 Sqrt[3]), {z, 1/6, 1/4}]

for i = 0, 1,...,14

and

Integrate[(4293331842800 z^i ArcTan[(2 - 8 z - 3 Sqrt[-1 + 16 z - 32z^2])/Sqrt[9 + 96 z^2 + 12 Sqrt[-1 + 16 z - 32 z^2] - 48 z (1 + Sqrt[-1 + 16 z - 32 z^2])]])/(729 Sqrt[2]), {z, 1/8 (2 - Sqrt[3]), 1/24 (6 - Sqrt[6])}]

also for i = 0, 1,...,14.

For these two latter sets of integrations, I needed to apply the denestSqrt function to the denominators of the arguments of the ArcTan function.

Summing all the results and performing a Together on the outcome, yielded (LeafCount=511)

(1/169075682574336)(194741930767680878400 [Pi] - 27558126743144112045 Sqrt2 [Pi] + 1635658501120 Sqrt[3] [Pi] - 1495234540800 Sqrt[6 (5 - 2 Sqrt[6])] [Pi] - 75180397824 Sqrt[2 (7 - 2 Sqrt[6])] [Pi] + 386620534272 Sqrt[3 (7 - 2 Sqrt[6])] [Pi] - 1495234540800 Sqrt[6 (5 + 2 Sqrt[6])] [Pi] + 75180397824 Sqrt[2 (7 + 2 Sqrt[6])] [Pi] + 386620534272 Sqrt[3 (7 + 2 Sqrt[6])] [Pi] - 337364402615426309760 ArcCsc[Sqrt[3]] + 55116262154719215450 Sqrt2 ArcCsc[Sqrt[3]] - 389483861535361756800 ArcSec[Sqrt[3]] + 55116262154719215450 Sqrt2 ArcSec[Sqrt[3]] + 2736717551038466595 Sqrt2 ArcSin[Sqrt2/3] - 52119458919935447040 ArcSin[1/Sqrt[3]] - 2736717551038466595 Sqrt2 ArcSin[1/Sqrt[3]] + 2736717551038466595 Sqrt2 ArcSin[1/3 Sqrt[1/3 (11 - 4 Sqrt[7])]] - 13002646487040 Sqrt2 ArcTan[1/Sqrt2] - 150360795648 Sqrt[2 (7 + 2 Sqrt[6])] ArcTan[(Sqrt2 + Sqrt[3])/Sqrt[7 - 2 Sqrt[6]]] - 773241068544 Sqrt[3 (7 + 2 Sqrt[6])] ArcTan[(Sqrt2 + Sqrt[3])/Sqrt[7 - 2 Sqrt[6]]] + 2990469081600 Sqrt[6 (5 + 2 Sqrt[6])] ArcTan[1/5 (2 + Sqrt[6]) Sqrt[7/2 + Sqrt[6]]] + 150360795648 Sqrt[2 (7 - 2 Sqrt[6])] ArcTan[(-Sqrt2 + Sqrt[3])/Sqrt[7 + 2 Sqrt[6]]] - 773241068544 Sqrt[3 (7 - 2 Sqrt[6])] ArcTan[(-Sqrt2 + Sqrt[3])/Sqrt[7 + 2 Sqrt[6]]] + 2990469081600 Sqrt[6 (5 - 2 Sqrt[6])] ArcTan[(-2 + Sqrt[6])/Sqrt[14 + 4 Sqrt[6]]])

Various manipulations including multiple uses of WolframAlpha led to the answer given at the outset. Near the completion, I used the substitution command

/. ArcTan[(1904 Sqrt[2])/5983] -> (5 Pi - 8 ArcCos[-1/3]) /. Log[1/6 + Sqrt[2/3] + 1/2 I Sqrt[11/9 - (4 Sqrt[2/3])/3]] -> I (2 Pi/3 - ArcCos[-1/3]) /. Log[1/6 (1 - 2 Sqrt[6] + I Sqrt[11 + 4 Sqrt[6]])] -> I (4 Pi/3 - ArcCos[-1/3])

A generalization of the question posed would replace the radius-squared constraint < 3/8 by a variable.

The radius-squared 1/3 case--of strong quantum-information-theoretic interest (eqs. (17), (18)) 1998paper--appears to be decidely simpler than the 3/8 case, the principal subject of the question above.

Then, the command

GenericCylindricalDecomposition[1 > x && x >= y && y >= z && z >= 1 -x - y - z >= 0 && x^2 + y^2 + (1 - x - y - z)^2 + z^2 < 1/3, {x,y,z}][[1]]

yields (LeafCount of 409 vs. 667 in the 3/8 case)

(1/4 < x < 1/3 && (1 - x)/3 < y < x && 1/2 (1 - x - y) < z < y) || (1/3 < x < 1/12 (3 + Sqrt[3]) && (((1 - x)/3 < y < (1 - x)/3 + 1/3 Sqrt[x - 2 x^2] && 1/2 (1 - x - y) < z < y) || ((1 - x)/3 + 1/3 Sqrt[x - 2 x^2] < y < x && 1/2 (1 - x - y) < z <  1/2 (1 - x - y) + Sqrt[-1 + 6 x - 9 x^2 + 6 y - 6 x y - 9 y^2]/(2 Sqrt[3])))) || (1/12 (3 + Sqrt[3]) < x < 1/2 && (((1 - x)/3 < y < (1 - x)/3 + 1/3 Sqrt[x - 2 x^2] && 1/2 (1 - x - y) < z < y) || ((1 - x)/3 + 1/3 Sqrt[x - 2 x^2] < y < (1 - x)/3 + 2/3 Sqrt[x - 2 x^2] && 1/2 (1 - x - y) < z < 1/2 (1 - x - y) + Sqrt[-1 + 6 x - 9 x^2 + 6 y - 6 x y - 9 y^2]/(2 Sqrt[3]))))

If we integrate the integrand

9081072000 (x - y)^2 (x - z)^2 (y - z)^2 (-1 + 2 x + y + z)^2 (-1 + x+ 2 y + z)^2 (-1 + x + y + 2 z)^2  .

subject to this result over $$z \in [0,1]$$, we get a two-part nontrivial result. Integrating the first part over $$y \in [0,1]$$, then $$x \in [0,1]$$, gives us

294529/2244806784

then, similarly, the second part yielding

(-294529 + 1122660 Sqrt[3] \[Pi])/2244806784

the two results summing to

(35 \[Pi])/(23328 Sqrt[3])  .

This approach seems more expeditious than directly attempting the three-fold integration using the result of the GenericCylindricalDecomposition command.

We will now investigate the case in which the radius squared is variable rather than fixed at 1/3 or 3/8, but these seems very challenging.

In pursuing this variable radius-squared case, it emerged that $$\frac{17}{50} \in [\frac{1}{3},\frac{3}{8}]$$ was a particular value of interest.

We were able to again obtain an exact value for the separability probability ($$\approx 0.00484591$$) associated with $$\frac{17}{50}$$.

This exact value--which we are attempting to simplify--took the form

(1/984150000000000000)(3036682376243712 [Pi] + 876531143656000 Sqrt[3] [Pi] - 3003692491800 Sqrt[6 (26 - 15 Sqrt[3])] [Pi] + 8485527480079051251360 Sqrt[2/(14 - 5 Sqrt[3])] [Pi] + 4899121574841721549728 Sqrt[6/(14 - 5 Sqrt[3])] [Pi] - 771412257557912814480 Sqrt[2 (14 + 5 Sqrt[3])] [Pi] - 445375074841312447104 Sqrt[6 (14 + 5 Sqrt[3])] [Pi] - 3003692491800 Sqrt[6 (26 + 15 Sqrt[3])] [Pi] + 26783430491537132034375 Sqrt2 ArcCsc[3] - 26783430491537132034375 Sqrt2 ArcCsc[9] + 26783430491537132034375 Sqrt2 ArcSin[1/9] - 26783430491537132034375 Sqrt2 ArcSin[1/3] - 6007384983600 Sqrt[6 (26 + 15 Sqrt[3])] ArcTan[(-1 - 3 Sqrt[3])/Sqrt[84 - 30 Sqrt[3]]] - 4555023564365568 ArcTan[1/4 (3 - Sqrt[3])] - 2629593430968000 Sqrt[3] ArcTan[1/4 (3 - Sqrt[3])] - 4555023564365568 ArcTan[1/4 (3 + Sqrt[3])] + 2629593430968000 Sqrt[3] ArcTan[1/4 (3 + Sqrt[3])] - 771412257557912814480 Sqrt[2 (14 - 5 Sqrt[3])] ArcTan[(3 Sqrt[2 (14 + 5 Sqrt[3])])/(9 - Sqrt[3])] + 445375074841312447104 Sqrt[6 (14 - 5 Sqrt[3])] ArcTan[(3 Sqrt[2 (14 + 5 Sqrt[3])])/(9 - Sqrt[3])] + 8485527480079051251360 Sqrt[2/(14 + 5 Sqrt[3])] ArcTan[(3 Sqrt[2 (14 + 5 Sqrt[3])])/(9 - Sqrt[3])] - 4899121574841721549728 Sqrt[6/(14 + 5 Sqrt[3])] ArcTan[(3 Sqrt[2 (14 + 5 Sqrt[3])])/(9 - Sqrt[3])] + 6007384983600 Sqrt[6 (26 - 15 Sqrt[3])] ArcTan[(-1 + 3 Sqrt[3])/Sqrt[84 + 30 Sqrt[3]]] - 2277511782182784 I Log[6] + 1314796715484000 I Sqrt[3] Log[6] + 2277511782182784 I Log[12] - 1314796715484000 I Sqrt[3] Log[12] - 2277511782182784 I Log[24] - 1314796715484000 I Sqrt[3] Log[24] + 2277511782182784 I Log[48] + 1314796715484000 I Sqrt[3] Log[48] + 2277511782182784 I Log[-9 (-2 + Sqrt[3])] - 1314796715484000 I Sqrt[3] Log[-9 (-2 + Sqrt[3])] - 4555023564365568 I Log[3 (-1 + Sqrt[3])] + 2629593430968000 I Sqrt[3] Log[3 (-1 + Sqrt[3])] - 4555023564365568 I Log[3 (1 + Sqrt[3])] - 2629593430968000 I Sqrt[3] Log[3 (1 + Sqrt[3])] + 2277511782182784 I Log[9 (2 + Sqrt[3])] + 1314796715484000 I Sqrt[3] Log[9 (2 + Sqrt[3])] + 8485527480079051251360 I Sqrt[2/(14 - 5 Sqrt[3])] Log[-1 + 3 Sqrt[3]] + 4899121574841721549728 I Sqrt[6/(14 - 5 Sqrt[3])] Log[-1 + 3 Sqrt[3]] - 771412257557912814480 I Sqrt[2 (14 + 5 Sqrt[3])] Log[-1 + 3 Sqrt[3]] - 445375074841312447104 I Sqrt[6 (14 + 5 Sqrt[3])] Log[-1 + 3 Sqrt[3]] + 8485527480079051251360 I Sqrt[2/(14 + 5 Sqrt[3])] Log[1 + 3 Sqrt[3]] - 4899121574841721549728 I Sqrt[6/(14 + 5 Sqrt[3])] Log[1 + 3 Sqrt[3]] - 4242763740039525625680 I Sqrt[2/(14 + 5 Sqrt[3])] Log[14 + 3 Sqrt[3]] + 2449560787420860774864 I Sqrt[6/(14 + 5 Sqrt[3])] Log[14 + 3 Sqrt[3]] - 8485527480079051251360 I Sqrt[2/(14 + 5 Sqrt[3])] Log[1 + 9 Sqrt[3]] + 4899121574841721549728 I Sqrt[6/(14 + 5 Sqrt[3])] Log[1 + 9 Sqrt[3]] + 4242763740039525625680 I Sqrt[2/(14 + 5 Sqrt[3])] Log[122 + 9 Sqrt[3]] - 2449560787420860774864 I Sqrt[6/(14 + 5 Sqrt[3])] Log[122 + 9 Sqrt[3]] - 8485527480079051251360 I Sqrt[2/(14 - 5 Sqrt[3])] Log[1/2 (-1 - 3 Sqrt[3] - I Sqrt[84 - 30 Sqrt[3]])] - 4899121574841721549728 I Sqrt[6/(14 - 5 Sqrt[3])] Log[1/2 (-1 - 3 Sqrt[3] - I Sqrt[84 - 30 Sqrt[3]])] + 771412257557912814480 I Sqrt[2 (14 + 5 Sqrt[3])] Log[1/2 (-1 - 3 Sqrt[3] - I Sqrt[84 - 30 Sqrt[3]])] + 445375074841312447104 I Sqrt[6 (14 + 5 Sqrt[3])] Log[1/2 (-1 - 3 Sqrt[3] - I Sqrt[84 - 30 Sqrt[3]])])

• For what it's worth, there is now ResourceFunction["RadicalDenest"] for denesting heuristics. Nov 11 '20 at 16:17
• Thanks, Daniel Lichtblau--will keep this in mind. Nov 11 '20 at 23:45