# Verify a conjectured formula for a modification of a 3D constrained integration successfully solved using Mathematica

I've acquired somewhat indirect--and not fully conclusive--evidence that the solution of a certain three-dimensional constrained integration takes the form

1/576 (8 - 6 Sqrt[2] + 15 Sqrt[2] \[Pi] - 48 Sqrt[2] ArcTan[Sqrt[2]])


$$\approx 0.00227243$$. My question here is whether this conjecture can be formally verified (needless to say, with the use of Mathematica).

The problem in question is a modification of the three-dimensional integration successfully and quite impressively solved by user JimB in the answer

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


$$\approx 0.00365826$$

The specific integration problem for which we seek to verify our candidate solution (following the notation of Original3Dproblem) is

Integrate[Boole[Subscript[λ, 1] > Subscript[λ, 2] && Subscript[λ, 2] > Subscript[λ, 3] &&Subscript[λ, 3] > 1 - Subscript[λ, 1] - Subscript[λ, 2] - Subscript[λ, 3] && Subscript[λ, 1] -Subscript[λ, 3] < 2 Sqrt[Subscript[λ, 2] (1 - Subscript[λ, 1] - Subscript[λ, 2] - Subscript[λ,3])]], {Subscript[λ, 3], 0, 1}, {Subscript[λ, 2], 0, 1}, {Subscript[λ, 1], 0, 1}]  .


The (unmodified) question Original3Dproblem was also posed in constrained form, but converted to an unconstrained form employing a transformation suggested by N. Tessore,

change = {Subscript[λ, 1] -> x/(1 + 2 x), Subscript[λ, 2] -> y/(1 + y) (1 + x)/(1 + 2 x),Subscript[λ, 3] -> z 1/(1 + y) (1 + x)/(1 + 2 x)},


which clearly remains applicable for the present (modified) question, leading to the transformed unconstrained problem at hand

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


also conjecturally yielding the formula given at the outset.

Although we have not (yet) been able to solve this problem directly, we have solved--using Mathematica--the associated 2D-integration for the boundary area of the convex set, modifying the inequality constraint

 Subscript[λ, 1] -Subscript[λ, 3] < 2 Sqrt[Subscript[λ, 2] (1 - Subscript[λ, 1] - Subscript[λ, 2] - Subscript[λ,3])]


to the equality constraint

 Subscript[λ, 1] -Subscript[λ, 3] == 2 Sqrt[Subscript[λ, 2] (1 - Subscript[λ, 1] - Subscript[λ, 2] - Subscript[λ,3])].


The solution to this 2D problem we did find to be

1/96 (8 - 6 Sqrt[2] + 15 Sqrt[2] \[Pi] - 48 Sqrt[2] ArcTan[Sqrt[2]])


$$\approx 0.013634585$$.

The key to obtaining our conjectured formula

1/576 (8 - 6 Sqrt[2] + 15 Sqrt[2] \[Pi] - 48 Sqrt[2] ArcTan[Sqrt[2]])


for which we seek verification here, is that we found the (area/volume) ratio of 0.013634585916219 to a numerical integration estimation (0.002272430980282073) of the solution to the 3D-problem to be 6.000000015193957, clearly pointing to an exact value of 6.

If the area/volume ratio is, in fact, 6, then this might serve as a useful clue in identifying the specific nature of the set in question, if it falls within known categories. (As simply an example, a three-dimensional ball of radius $$\frac{1}{2}$$ has such a ratio.)

The modification pursued here consists in replacing the (Hilbert-Schmidt [eq. (15.35)] GeometryQuantumStates) integrand in Original3Dproblem

9081072000 (Subscript[λ, 1] - Subscript[λ, 2])^2 (Subscript[λ, 1] - Subscript[λ, 3])^2 (Subscript[λ, 2] - Subscript[λ, 3])^2 (-1 + 2 Subscript[λ, 1] + Subscript[λ, 2] + Subscript[λ, 3])^2 (-1 + Subscript[λ, 1] + 2 Subscript[λ, 2] + Subscript[λ, 3])^2 (-1 + Subscript[λ, 1] + Subscript[λ, 2] + 2 Subscript[λ, 3])^2


by simply 1.

The motivation behind this modification was that rather than considering the problem as that of the four ordered eigenvalues of a (Hermitian, nonnegative-definite $$4 \times 4$$, trace 1) "two-qubit density matrix" in the 15-dimensional setting for such such matrices, we simply now focus on the 3-dimensional convex set $$(\lambda_1, \lambda_2, \lambda_3, 1-\lambda_1-\lambda_2-\lambda_3)$$ of "ordered spectra of absolutely separable two-qubit density matrices".

We are interested in this problem because its solution would yield the Euclidean volume of that indicated convex set for which we aspire JohnEllipsoidProblem to find the "John ellipsoids" of minimal and maximal volumes circumscribing and inscribing it.

• Is this a question about Mathematica? This seems like it might be more appropriate for Mathematics, maybe Physics. It could also use some editing to make clear what exactly you are asking for. Nov 16, 2020 at 21:36
• N.J.Evans--thanks! I've done a lot or reorganizing to try to bring the Mathematica aspects to the fore, and putting the motivation/background later in the question statement. Nov 17, 2020 at 4:05

For your first question $$\frac{1}{2} \cos ^{-1}\left(\frac{1}{3}\right)-\frac{\pi }{8}$$ is equivalent to $$\csc ^{-1}\left(\sqrt{6 \left(\sqrt{2}+2\right)}\right)$$ so the equation can be simplified to

1/288 (4 - 3 Sqrt[2] - 6 Sqrt[2] ArcCsc[3] + 12 Sqrt[2] ArcCsc[Sqrt[6 (2 + Sqrt[2])]]) /.
ArcCsc[Sqrt[6 (2 + Sqrt[2])]] -> -(π/8) + 1/2 ArcCos[1/3] /.
ArcCsc[3] -> π/2 - ArcCos[1/3] // Expand // Together

(*  1/576 (8 - 6 Sqrt[2] - 9 Sqrt[2] π + 24 Sqrt[2] ArcCos[1/3]) *)


as in previous questions you seemed to desire the term ArcCos[1/3} to be included.

The next part is to use Mathematica to end up with that result.

Taking the Boole part of the formula one can end up with 5 integrations to be performed:

Reduce[Subscript[λ, 1] > Subscript[λ, 2] && Subscript[λ, 2] > Subscript[λ, 3] &&
Subscript[λ, 3] > 1 - Subscript[λ, 1] - Subscript[λ, 2] - Subscript[λ, 3] &&
Subscript[λ, 1] - Subscript[λ, 3] < 2 Sqrt[Subscript[λ, 2] (1 - Subscript[λ, 1] - Subscript[λ, 2] - Subscript[λ, 3])]]


And after some manipulations of that result one ends up with 5 regions to integrate over:

Next one performs the individual integrations:

Integration 1

i1 = Integrate[1, {λ1, 1/4, 1/3}, {λ2, 1/3 (1 - λ1), λ1}, {λ3, 1/2 (1 - λ1 - λ2), λ2}]
(* 1/3888 *)


Integration 2:

i2 = Integrate[1, {λ1, 1/3, 1/8 (2 + Sqrt[2])},
{λ2, 1/3 (1 - λ1), 1/9 (2 - λ1) + 2/9 Sqrt[1 - λ1 - 2 λ1^2]},
{λ3, 1/2 (1 - λ1 - λ2), λ2}]
(* (998-447 Sqrt[2]+972 Sqrt[2] ArcSin[1/54 (20-23 Sqrt[2])])/62208 *)


The 3rd integration results in some complicated square root functions and the denestSqrt function from @CarlWoll is used.

denestSqrt[e_, domain_, x_] := Replace[y /. Solve[Simplify[Reduce[Reduce[y == e && domain, x], y, Reals], domain], y], {{r_} :> r, _ -> e}]

i3a = Integrate[1, {λ2, 1/9 (2 - λ1) + 2/9 Sqrt[1 - λ1 - 2 λ1^2], λ1},
{λ3, 1/2 (1 - λ1 - λ2), λ1 - 2 λ2 + 2 Sqrt[λ2 - 2 λ1 λ2]},
Assumptions -> {1/3 < λ1 <= 1/8 (2 + Sqrt[2])}] // Expand;
i3a = i3a /. Sqrt[1 - 2 λ1] Sqrt[2 - λ1 + 2 Sqrt[1 - λ1 - 2 λ1^2]] ->
denestSqrt[Sqrt[(1 - 2 λ1) (2 - λ1 + 2 Sqrt[1 - λ1 - 2 λ1^2])], 1/3 < λ1 <= 1/8 (2 + Sqrt[2]), λ1] /.
Sqrt[1 - 2 λ1] λ1 Sqrt[2 - λ1 + 2 Sqrt[1 - λ1 - 2 λ1^2]] ->
denestSqrt[λ1 Sqrt[(1 - 2 λ1) (2 - λ1 + 2 Sqrt[1 - λ1 - 2 λ1^2])],
1/3 < λ1 <= 1/8 (2 + Sqrt[2]), λ1] /.
Sqrt[1 - 2 λ1] Sqrt[1 - λ1 - 2 λ1^2] Sqrt[2 - λ1 + 2 Sqrt[1 - λ1 - 2 λ1^2]] ->
denestSqrt[Sqrt[(1 - 2 λ1) (1 - λ1 - 2 λ1^2) (2 - λ1 + 2 Sqrt[1 - λ1 - 2 λ1^2])],
1/3 < λ1 <= 1/8 (2 + Sqrt[2]), λ1] // Expand;
i3a1 = Integrate[-(1/81), {λ1, 1/3, 1/8 (2 + Sqrt[2])}];
i3a2 = Integrate[-((50 λ1)/81), {λ1, 1/3, 1/8 (2 + Sqrt[2])}];
i3a3 = Integrate[4/3 Sqrt[1 - 2 λ1] λ1^(3/2), {λ1, 1/3, 1/8 (2 + Sqrt[2])}] // ToRadicals;
i3a4 = Integrate[(77 λ1^2)/81, {λ1, 1/3, 1/8 (2 + Sqrt[2])}];
i3a5 = Integrate[-(1/81) Sqrt[1 - λ1 - 2 λ1^2], {λ1, 1/3, 1/8 (2 + Sqrt[2])}];
i3a6 = Integrate[-(10/81) λ1 Sqrt[1 - λ1 - 2 λ1^2], {λ1, 1/3, 1/8 (2 + Sqrt[2])}];
i3 = i3a1 + i3a2 + i3a3 + i3a4 + i3a5 + i3a6 // Expand
(* -(329/31104)+133/(31104 Sqrt[2])-ArcSin[1/54 (20-23 Sqrt[2])]/(96 Sqrt[2])+ArcSin[1/2 Sqrt[1/3 (2-Sqrt[2])]]/(24 Sqrt[2]) *)


Integration 4

i4 = Integrate[1, {λ1, 1/8 (2 + Sqrt[2]), 1/2},
{λ2, 1/3 (1 - λ1), 1/9 (2 - λ1) + 2/9 Sqrt[1 - λ1 - 2 λ1^2]},
{λ3, 1/2 (1 - λ1 - λ2), λ2}]
(* (-2+149 Sqrt[2]-324 Sqrt[2] ArcCos[1/6 (4+Sqrt[2])])/20736 *)


Integration 5:

i5a = Integrate[1, {λ2, 1/9 (2 - λ1) + 2/9 Sqrt[1 - λ1 - 2 λ1^2],
1/9 (5 - 7 λ1) + 4/9 Sqrt[1 - λ1 - 2 λ1^2]},
{λ3, 1/2 (1 - λ1 - λ2), λ1 - 2 λ2 + 2 Sqrt[λ2 - 2 λ1 λ2]},
Assumptions -> {1/8 (2 + Sqrt[2]) < λ1 < 1/2}] // Expand;

i5a = i5a /. Sqrt[1 - 2 λ1] Sqrt[2 - λ1 + 2 Sqrt[1 - λ1 - 2 λ1^2]] ->
denestSqrt[Sqrt[(1 - 2 λ1) (2 - λ1 + 2 Sqrt[1 - λ1 - 2 λ1^2])], 1/8 (2 + Sqrt[2]) < λ1 < 1/2, λ1] /.
Sqrt[1 - 2 λ1] λ1 Sqrt[2 - λ1 + 2 Sqrt[1 - λ1 - 2 λ1^2]] ->
λ1 denestSqrt[Sqrt[(1 - 2 λ1) (2 - λ1 + 2 Sqrt[1 - λ1 - 2 λ1^2])], 1/8 (2 + Sqrt[2]) < λ1 < 1/2, λ1] /.
Sqrt[1 - 2 λ1] Sqrt[1 - λ1 - 2 λ1^2] Sqrt[2 - λ1 + 2 Sqrt[1 - λ1 - 2 λ1^2]] ->
denestSqrt[Sqrt[(1 - 2 λ1) (1 - λ1 - 2 λ1^2) (2 - λ1 + 2 Sqrt[1 - λ1 - 2 λ1^2])], 1/8 (2 + Sqrt[2]) < λ1 < 1/2, λ1] /.
Sqrt[1 - 2 λ1] Sqrt[5 - 7 λ1 + 4 Sqrt[1 - λ1 - 2 λ1^2]] ->
denestSqrt[Sqrt[(1 - 2 λ1) (5 - 7 λ1 + 4 Sqrt[1 - λ1 - 2 λ1^2])], 1/8 (2 + Sqrt[2]) < λ1 < 1/2, λ1] /.
Sqrt[1 - 2 λ1] λ1 Sqrt[5 - 7 λ1 + 4 Sqrt[1 - λ1 - 2 λ1^2]] ->
λ1 denestSqrt[Sqrt[(1 - 2 λ1) (5 - 7 λ1 + 4 Sqrt[1 - λ1 - 2 λ1^2])], 1/8 (2 + Sqrt[2]) < λ1 < 1/2, λ1] /.
Sqrt[1 - 2 λ1] Sqrt[1 - λ1 - 2 λ1^2] Sqrt[5 - 7 λ1 + 4 Sqrt[1 - λ1 - 2 λ1^2]] ->
denestSqrt[Sqrt[(1 - 2 λ1) (1 - λ1 - 2 λ1^2) (5 - 7 λ1 + 4 Sqrt[1 - λ1 - 2 λ1^2])], 1/8 (2 + Sqrt[2]) < λ1 < 1/2, λ1] // Expand
(* 7/324+(2 λ1)/81-(11 λ1^2)/81+1/27 Sqrt[1-λ1-2 λ1^2]-2/27 λ1 Sqrt[1-λ1-2 λ1^2] *)
i5 = Integrate[i5a, {λ1, 1/8 (2 + Sqrt[2]), 1/2}]
(* (514-781 Sqrt[2]+972 Sqrt[2] ArcCos[1/6 (4+Sqrt[2])])/62208 *)


result = i1 + i2 + i3 + i4 + i5 // FullSimplify
1/288 (4 - 3 Sqrt[2] + 6 Sqrt[2] ArcCsc[Sqrt[6 (2 + Sqrt[2])]] + 3 Sqrt[2] ArcSin[1/54 (20 - 23 Sqrt[2])])


This can be further simplified to

result /. ArcSin[1/54 (20 - 23 Sqrt[2])] -> -((5 \[Pi])/4) + 3 ArcCos[1/3] /.
ArcCsc[Sqrt[6 (2 + Sqrt[2])]] -> -(\[Pi]/8) + 1/2 ArcCos[1/3] // Expand // Together
(* 1/576 (8 - 6 Sqrt[2] - 9 Sqrt[2] \[Pi] + 24 Sqrt[2] ArcCos[1/3]) *)

• Thanks, JimB! I probably should have--for full consistency--employed this simplification, but I used a slight variation with ArcTan that I introduced in a revision of the original question. Nov 17, 2020 at 4:02
• Fabulous, JimB, as usual! Not too easy once again--despite the severe simplification of the integrand. Regarding the slight difference in final formulas (in terms of inverse trigonometric formulas), when I use WolframAlpha for exact results based on high-precision numerics, it seems to always suggest alternatives based on use of the "magic angle" or "tetrahedral angle". This time I happened to pick the magic one--but probably for overall consistency, I should have picked the alternative. Nov 17, 2020 at 13:13
• Of course, one has better things to do--but I think it would be an interesting exercise to try to solve the problem in its original, constrained form. This must have been the way I pursued the original/unmodified challenging problem in 2009--the Tessore transformation not being available then. (The last answer to mathematica.stackexchange.com/questions/230127/… did pursue the constrained problem.) Nov 17, 2020 at 13:22