# Evaluate a certain three-dimensional constrained integral

The result of the three-dimensional integration

Integrate[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 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}],


that is,

for the two-qubit Hilbert-Schmidt absolute separability probability apparently can be expressed as

$$$$\label{HSabs} \frac{29902415923}{497664}+\frac{-3217542976+5120883075 \pi -16386825840 \tan ^{-1}\left(\sqrt{2}\right)}{32768 \sqrt{2}} =$$$$ $$$$\frac{32(29902415923 - 24433216974 \sqrt{2})+248874917445 \sqrt{2}(5 \pi - 16 \tan ^{-1}\left(\sqrt{2}\right))}{2^{16} \cdot 3^5} \approx 0.00365826$$$$

QuantumComputingStackExchangeQuestion

Can this be explicitly confirmed using Mathematica?

Through use of the transformation,

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)};


Nicolas Tessore has now reported to me that he was able to convert the 3D integral into an unconstrained one of the form,

Integrate[integrand2, {z, 1/2, 1}, {y, z, 2 + 2 Sqrt[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}],


where

integrand2 = (9081072000 (1 + x)^8 (x - y)^2 (1 - 2 z)^2 (y - z)^2 (-1 + y + z)^2 (z + x (-1 - y + z))^2 (-1 + z + x (y + z))^2)/((1 + 2 x)^16 (1 + y)^15)  .


Let me indicate here that the indicated result was obtained in the 2009paper

2009paper

I obtained this result (eq. (34) there), but the now-requested step-by-step process was not detailed. Comments of present interest there were that [C]opious use was made of trigonometric identities involving the tetrahedral dihedral angle $$\phi=\cos ^{-1}\left(\frac{1}{3}\right)$$, assisted by V. Jovovic" and that use was made of the Sloane website sequence A025172--"Let phi = arccos(1/3), the dihedral angle of the regular tetrahedron. Then cos(n*phi) = a(n)/3^n". (This sequence is "[u]sed when showing that the regular simplex is not "scisssors-dissectable" to a cube, thus answering Hilbert's third problem.".)

A025172sequence

These comments led me to consult my email archives.

On April 21, 2008 I wrote to Vladeta Jovovic (and also Wouter Meeussen and Neil Sloane) the following:

"Dear Drs. Jovocic/Meeussen/Sloane:

I been doing some analyses in which I've been able to simplify several trigonometric terms using the relation

cos(n*phi) =a(n)/3^n where phi=ArcCos(1/3)

reported in Sloane's Superseeker sequence

A025172.

I have two further terms ArcTan[(1/329 )(729 - 460 Sqrt[2])] and ArcTan[(1/7) (9 + 4 Sqrt[2])], which also clearly pertain, since 329/729 =a(6)/3^6 and -7/9 =a(2)/3^2.

But I don't see how to manipulate them to reexpress/simplify them in terms of phi, which I presume is possible/natural.

Perhaps you have some insights in this matter?

Sincerely,

Paul B. Slater

P. S. I also have the term

ArcTan[(1/7) (-3 + Sqrt[2])]


which perhaps also has some simpler form."

"for n from 0 to 10 do q:=tan(-n*phi):print(expand(q));od:

                                0

1/2
-2 2

1/2
4 2
------
7

1/2
10 2
- -------
23

1/2
56 2
-------
17

1/2
22 2
-------
241

1/2
460 2
- --------
329

1/2
1118 2
---------
1511

1/2
1904 2
- ---------
5983

1/2
13870 2
----------
1633

1/2
10604 2
----------
57113

V.



and

phi=ArcCos(1/3)

ArcTan[(1/329 )(729 - 460 Sqrt[2])]


5Pi/4 - 3phi

ArcTan[(1/7) (9 + 4 Sqrt[2])],


3*Pi/4 - phi.

Within the next week, V. Jovovic also wrote:

ArcTan[(1/7) (-9 + 4 Sqrt[2])]


Pi/4-phi

ArcTan[(1/7) (-3 + Sqrt[2])]


Pi/8-phi/2

ArcSin[(1/6) (4 + Sqrt[2])]
= 3*Pi/4 - phi


and

ArcCsc[3/17 Sqrt[52 + 14 Sqrt[2]]]


5*Pi/8-phi

ArcTan[7/(3 + Sqrt[2])]


Pi/8+phi/2

ArcTan[1/(3 + Sqrt[2])]

• Pi/8+phi/2
ArcCsc[Sqrt[6 (2 + Sqrt[2])]]


5*Pi/8-phi

Although this 2008 email correspondence was clearly central to the obtaining of the indicated formula (for which a Mathematica demonstration is requested), it is presently not clear to me in what manner the results discussed there were obtained and further employed. (Also, apparently this 2008 correspondence was carried on after(!) I had been able to perform the desired 3D integration, and had a result for which some simplification--using the Jovovic transformations--was possible.)

• If you substitude Integrate by NIntegrate you'll get the numerical verification ...=0.00365826  Sep 15 '20 at 8:01
• @UlrichNeumann It is not explicitly confirmed, so try to use Integrate. Sep 15 '20 at 17:15
• @AlexTrounev I tried Integrate but it doesn't evaluate. Additional I checked the integration range (defined by Block) which is composed by three planes. Perhaps splitting Integrate helps solving symbolically. Sep 15 '20 at 17:55
• The command LeafCount[Integrate[the integrand under consideration,{Subscript[\[Lambda], 1], 0, 1}, Assumptions -> Subscript[\[Lambda], 2] >= 0 && Subscript[\[Lambda], 2] <= 1 && Subscript[\[Lambda], 3] >= 0 && Subscript[\[Lambda], 3] <= 1]] performs 5177. To much for a next integration. Sep 16 '20 at 15:33
• Not sure it helps much but 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 in the updated Integrate[integrand2...] statement can be simplified to (4 y + z - 3 y z - z^2 + 2 (1 + y) Sqrt[y - y z])/(-1 + y + z)^2.
– JimB
Sep 16 '20 at 23:33

This is now closer to an answer in that I attempted to follow @UlrichNeumann 's good suggestion about splitting the integration into parts.

I changed the subscripted variables to x, y, and z to lighten the text load. The constraints in the Boole function can be written as

1 > x > y > z > 0 && z > 1 - x - y - z && x - z < 2 Sqrt[y (1 - x - y - z)]


If Reduce is used on this

Reduce[1 > x > y > z > 0 && z > 1 - x - y - z && x - z < 2 Sqrt[y (1 - x - y - z)]]


one can see 5 mutually exclusive integrations to be performed:

c1 = 1/8 (2 - Sqrt[2]) < z <= 1/6 && 1 - z - 2 Sqrt[z - 2 z^2] < y <= 1/2 (1 - 2 z) &&
1 - y - 2 z < x < -2 y + z + 2 Sqrt[y - 2 y z];
c2 = 1/8 (2 - Sqrt[2]) < z <= 1/6 && 1/2 (1 - 2 z) < y < (2 - z)/9 + 2/9 Sqrt[1 - z - 2 z^2] &&
y < x < -2 y + z + 2 Sqrt[y - 2 y z];
c3 = 1/6 < z <= 1/4 && z < y <= 1/2 (1 - 2 z) && 1 - y - 2 z < x < -2 y + z + 2 Sqrt[y - 2 y z];
c4 = 1/6 < z <= 1/4 && 1/2 (1 - 2 z) < y < (2 - z)/9 + 2/9 Sqrt[1 - z - 2 z^2] &&
y < x < -2 y + z + 2 Sqrt[y - 2 y z];
c5 = 1/4 < z < 1/3 && z < y < (2 - z)/9 + 2/9 Sqrt[1 - z - 2 z^2] &&
y < x < -2 y + z + 2 Sqrt[y - 2 y z];


The integrations over x and y can be written as follows:

(* Define the integrand *)
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;

(* Initial integrations *)
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}];(* Around 6 minutes *)
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}]; (* Around 4 minutes *)
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}]; (* Around 3 minutes *)
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]},
Assumptions -> {1/4 < z < 1/3}]; (* Around 5 minutes *)


Now some simplifications are needed for the integrals that Mathematica and Rubi don't handle directly:

(* Simplification rules *)
(* These can be shown to be true (or at least reasonable) by looking at Series[...,{z,0,10}] *)
(* These are the terms that Mathematica and Rubi can't handle directly *)
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]};


Combine into 3 integrands with the same range of values for z

integrandA = Expand[Expand[a1 + a2] /. rules] /. Sqrt[(1 - 2 z) z^k_] -> z^((k - 1)/2) Sqrt[(1 - 2 z) z];
integrandB = Expand[Expand[a3 + a4] /. rules] /. Sqrt[(1 - 2 z) z^k_] -> z^((k - 1)/2) Sqrt[(1 - 2 z) z];
integrandC = Expand[Expand[a5] /. rules] /. Sqrt[(1 - 2 z) z^k_] -> z^((k - 1)/2) Sqrt[(1 - 2 z) z];



All 3 integrands have a common structure and shared coefficients that can be written as follows:

$$\text{Integral of integrandA}=\sum_{i=0}^{14} c_{A,i}\int\limits_{\frac{2-\sqrt{2}}{8}}^{1/6}z^i dz+ \sum_{i=5}^{13}d_i\int\limits_{\frac{2-\sqrt{2}}{8}}^{1/6}z^i \sqrt{(1-2z)z}dz+ \sum_{i=5}^{13}e_i\int\limits_{\frac{2-\sqrt{2}}{8}}^{1/6}z^i \sqrt{1-z-2z^2}dz$$

$$\text{Integral of integrandB}=\sum_{i=0}^{14} c_{B,i}\int\limits_{\frac{2-\sqrt{2}}{8}}^{1/6}z^i dz- \sum_{i=5}^{13}d_i\int\limits_{\frac{2-\sqrt{2}}{8}}^{1/6}z^i \sqrt{(1-2z)z}dz+ \sum_{i=5}^{13}e_i\int\limits_{\frac{2-\sqrt{2}}{8}}^{1/6}z^i \sqrt{1-z-2z^2}dz$$

$$\text{Integral of integrandC}=\sum_{i=0}^{14} c_{C,i}\int\limits_{\frac{2-\sqrt{2}}{8}}^{1/6}z^i dz- \sum_{i=5}^{13}d_i\int\limits_{\frac{2-\sqrt{2}}{8}}^{1/6}z^i \sqrt{(1-2z)z}dz+ \sum_{i=5}^{13}e_i\int\limits_{\frac{2-\sqrt{2}}{8}}^{1/6}z^i \sqrt{1-z-2z^2}dz$$

So because some of the coefficients are shared one can simplify the total a bit by performing integrations that go from {z,1/6, 1/3} and {z, 1/8 (2 - Sqrt[2]), 1/3}. The reason for doing so is that the result is in a much simpler form. The following is the final integral:

$$\begin{multline} \sum_{i=0}^{14} c_{A1,i}\int\limits_{\frac{2-\sqrt{2}}{8}}^{1/6}z^i dz + \sum_{i=0}^{14} c_{B,i}\int\limits_{1/6}^{1/4}z^i dz + \sum_{i=0}^{14} c_{C,i}\int\limits_{1/4}^{1/3}z^i dz + \\ \sum_{i=5}^{13}d_i\int\limits_{\frac{2-\sqrt{2}}{8}}^{1/6}z^i \sqrt{(1-2z)z}dz- \sum_{i=5}^{13}d_i\int\limits_{1/6}^{1/3}z^i \sqrt{(1-2z)z}dz+ \sum_{i=5}^{13}e_i\int\limits_{\frac{2-\sqrt{2}}{8}}^{1/3}z^i \sqrt{1-z-2z^2}dz \end{multline}$$

In terms of Mathematica (which might show my lack of skill in extracting the appropriate pieces):

(* Now collect integrand terms to allow for combining ranges of z *)
i1 = integrandA /. Sqrt[1 - z - 2 z^2] -> 0 /. Sqrt[(1 - 2 z) z] -> 0;
i2 = integrandB /. Sqrt[1 - z - 2 z^2] -> 0 /. Sqrt[(1 - 2 z) z] -> 0;
i3 = integrandC /. Sqrt[1 - z - 2 z^2] -> 0 /. Sqrt[(1 - 2 z) z] -> 0;
i4 = integrandA - (integrandA /. Sqrt[(1 - 2 z) z] -> 0);
i5 = integrandA - (integrandA /. Sqrt[1 - z - 2 z^2] -> 0);
s1 = Integrate[i1, {z, 1/8 (2 - Sqrt[2]), 1/6}] // Expand // FullSimplify;
s2 = Integrate[i2, {z, 1/6, 1/4}] // Expand // FullSimplify;
s3 = Integrate[i3, {z, 1/4, 1/3}] // Expand // FullSimplify;
s4 = Integrate[i4, {z, 1/8 (2 - Sqrt[2]), 1/6}] // Expand // FullSimplify;
s5 = -Integrate[i4, {z, 1/6, 1/3}] // Expand // FullSimplify;
s6 = Integrate[i5, {z, 1/8 (2 - Sqrt[2]), 1/3}] // Expand // FullSimplify;
total = s1 + s2 + s3 + s4 + s5 + s6 // FullSimplify

(* (746624752335 Sqrt[2] π - 4 (-478438654768 + 390931471584 Sqrt[2] +
497749834890 Sqrt[2] ArcCsc[3] +
497749834890 Sqrt[2] ArcSec[Sqrt[3]] -
248874917445 Sqrt[2] ArcSin[1/6 (4 - Sqrt[2])]))/31850496 *)

N[s1 + s2 + s3 + s4 + s5 + s6, 50]
(* 0.0036582630543034854603976004088368426270721405774277 *)


A bit more simplification for the trigonometric terms:

ArcCsc[3] = π/2-ϕ where ϕ=ArcCos[1/3], ArcSec[Sqrt[3]]=π/2 - ϕ/2, and ArcSin[1/6 (4 - Sqrt[2])] = -π/4 + ϕ. That simplifies the result to

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


or

$$\frac{29902415923}{497664}-\frac{50274109}{512 \sqrt{2}}-\frac{3072529845 \pi }{32768 \sqrt{2}}+\frac{1024176615 \cos ^{-1}\left(\frac{1}{3}\right)}{4096 \sqrt{2}}$$

• Wow--that update really looks impressive! The $k=2$ formula with which it ends contains the Sqrt[2]'s that appears in the presumed full answer given in my question, and if one adds ArcSin[1/3] to ArcCsc[Sqrt[3]], one gets the ArcTan[Sqrt[2]] appearing there too. In the question I gave details of my 2008 correspondence with V. Jovovic focused on the Sloane Superseeker sequence A025172, pertaining to the clearly related phi = ArcCos[1/3], which was used repeatedly in the lengthy simplification process that led to the presumed answer given in the question. So, this all looks very right-headed! Sep 20 '20 at 12:16
• I've now solved all of the 5 integrals and the result gives the same number as the numerical approximation. I'll include that in an update tonight or tomorrow. However, there are ArcSin terms that I haven't yet been able to simplify. The result I get so far is....
– JimB
Sep 21 '20 at 23:04
• 29902415923/497664 - 50274109/(512 Sqrt[2]) + (48882576431592402695 ArcSin[5/9])/(194055791689728 Sqrt[2]) + (1024176615 ArcSin[7/9])/(16384 Sqrt[2]) - ( 1024176615 ArcSin[Sqrt[2/3]])/(8192 Sqrt[2]) + (1024176615 ArcSin[1/Sqrt[3]])/(8192 Sqrt[2]) + (1024176615 ArcSin[1/2 Sqrt[1/3 (2 - Sqrt[2])]])/(8192 Sqrt[2]) - ( 1861973732913175 Sqrt[2] ArcSin[1/6 (4 - Sqrt[2])])/11844225567 - ( 48882576431592402695 ArcSin[1/54 (20 + 5 Sqrt[2] + 4 Sqrt[7] - 8 Sqrt[14])])/(194055791689728 Sqrt[2]).
– JimB
Sep 21 '20 at 23:06
• Fabulous,JimB! I'm still mired in attacking the original, constrained (pre-Tessore contribution) integration problem. Having to develop rules for converting integration limits involving square roots of square roots into limits involving sums of single square roots--so the integrations are performable. Some (but not all) of those rules were given in the answers to math.stackexchange.com/questions/3375075/… . So, your current result involves only Sqrt[2] and ArcSin's--very neat! Sep 22 '20 at 12:00
• I see you also have Sqrt[7] and Sqrt[14]. Commenting further on my (lagging) analysis, I have a integration limit transformation rule of the form Sqrt[(1 - 2 x) (2 - x + 2 Sqrt[1 - x - 2 x^2])] -> 1 - 2 x + Sqrt[1 - x - 2 x^2], which Mathematica implements finely, but I also need Sqrt[(-1 + 2 x) (2 - x + 2 Sqrt[1 - x - 2 x^2])] -> 1 - 2 x + Sqrt[1 - x - 2 x^2],--note changes of sign in first Sqrt--which Mathematica immediately converts to the previous rule, so the second rule doesn't get enforced. So, I need to suppress the immediate conversion. Sep 22 '20 at 12:16

Well (by no means an actual answer), here's an interesting first step to what appears to be a very daunting (but apparently eventually solvable) problem. (Also, stemming from current social restrictions, my personal/home computer resources are on the minimal side.)

If we perform

T = CoefficientList[Numerator[integrand2], x]

where, as indicated in the question,

    integrand2 = (9081072000 (1 + x)^8 (x - y)^2 (1 - 2 z)^2 (y - z)^2 (-1 + y + z)^2 (z + x (-1 - y + z))^2 (-1 + z + x (y + z))^2)/((1 + 2 x)^16 (1 + y)^15)


is the unconstrained integrand provided by N. Tessore, we obtain

Dimensions[T][[1]] = 15 .


Now, using the simplified integration limits over $$x$$

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


provided by JimB in his comment to the question, the command (consulting the full structure--specifically, the denominator--of integrand2)

W = Integrate[x^i/(1 + 2 x)^16, {x, y, (4 y + z - 3 y z - z^2 +
2 (1 + y) Sqrt[y - y z])/(-1 + y + z)^2}, Assumptions -> 1/2 < z < 1 && z < y < 2 + 2 Sqrt[1 - z] - z && i >= 0]


interestingly yields

R  =    (1/1307674368000)(( 638512875 y^(-15 + i) Hypergeometric2F1[16, 15 - i, 16 - i, -(1/(2 y))])/(32 (15 - i)) + (1/(32 (15 - i))) 638512875 E^(-I i \[Pi]) (-1 + y + z)^(30 - 2 i) (-4 y - 2 y^(3/2) Sqrt[1 - z] - 2 Sqrt[y (1 - z)] - z + 3 y z + z^2)^(-15 + i)Hypergeometric2F1[16, 15 - i, 16 - i, (-1 + y + z)^2/(2 (-4 y - 2 y^(3/2) Sqrt[1 - z] - 2 Sqrt[y (1 - z)] - z + 3 y z + z^2))])


Using Simplify (I'm not sure if FullSimplify would be productive--but only tried briefly), this becomes

 S =   (1/(65536 (15 - i)))(y^(-15 + i) Hypergeometric2F1[16, 15 - i, 16 - i, -(1/(2 y))] + E^(-I i \[Pi]) (-1 + y + z)^(30 - 2 i) (-2 y^(3/2) Sqrt[1 - z] - z + z^2 + y (-4 + 3 z) - 2 Sqrt[y - y z])^(-15 + i)Hypergeometric2F1[16, 15 - i, 16 - i, -((-1 + y + z)^2/(y (8 - 6 z) + 4 y^(3/2) Sqrt[1 - z] + 2 z - 2 z^2 + 4 Sqrt[y - y z]))])


For $$i=0,\ldots,14$$, we obtain, it appears, rational functions--which need to be summed--of $$y$$ and $$z$$.

So, conceptually at least,

Sum[T[[i]] S/(1+y)^{15}, {i, 0, 14}]


yields the result of the outer of the required three integrations. I just did the indicated sum [LeafCount of 131971], followed by a Together [LeafCount 61470], with denominator

40040 (1 + 2 y)^15 (1 + 6 y + y^2 + 4 y^(3/2) Sqrt[1 - z] + 4 Sqrt[-y (-1 + z)] - 4 y z - z^2)^15


Simplify on the numerator gives LeafCount of 28806. FullSimplify on this result only gives LeafCount of 28770.

Extending the above approach, I just tried

WW = Integrate[x^i y^j/((1 + y)^(15) (1 + 2 x)^(16)), {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}, Assumptions -> 1/2 < z < 1 && i >= 0 && j >= 0]


to see if the two--as opposed to just one--outer integrations could be similarly performed. But unsuccessful (returned unevaluated).

Could also try the last twofold integration for specific values of i and j, $$i=0,\ldots,14$$, $$j=0,\ldots,10$$.

I just had a bizarre occurrence. I started a new kernel, entered integrand2, and

 T   = CoefficientList[Numerator[integrand2], x]


as previously, and for some unknown reason, I got a PolynomialReduce (a command I'd never used, and really don't understand) result of

{{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, 9081072000 y^4 z^2-18162144000 y^5 z^2 + 9081072000 y^6 z^2 - 18162144000 y^3 z^3 -36324288000 y^4 z^3 + 108972864000 y^5 z^3 - 54486432000 y^6 z^3 +9081072000 y^2 z^4 + 127135008000 y^3 z^4 - 9081072000 y^4 z^4 -236107872000 y^5 z^4 + 118053936000 y^6 z^4 - 72648576000 y^2 z^5 -345080736000 y^3 z^5 + 236107872000 y^4 z^5 + 217945728000 y^5 z^5 -108972864000 y^6 z^5 + 236107872000 y^2 z^6 + 454053600000 y^3 z^6 -417729312000 y^4 z^6 - 72648576000 y^5 z^6 + 36324288000 y^6 z^6 -399567168000 y^2 z^7 - 290594304000 y^3 z^7 + 290594304000 y^4 z^7 + 372323952000 y^2 z^8 + 72648576000 y^3 z^8 - 72648576000 y^4 z^8 - 181621440000 y^2 z^9 + 36324288000 y^2 z^10}


????????

The answer of JimB relied upon the transformation obtained by N. Tessore by which the originally constrained integration problem is convertable into an unconstrained one.

Since the original analysis reported in 2009paper proceeded with the constrained problem (the Tessore transformation only being reported recently), it remained a matter of some curiosity--at least to me--as to how one might proceed remaining within the constrained framework.

So proceeding, I issued the command (using the transformation {Subscript[\[Lambda], 1] -> x, Subscript[\[Lambda], 2] -> y, Subscript[\[Lambda], 3] -> z} on the original variables),

GenericCylindricalDecomposition[x > y && y > z && z > 1 - x - y - z && x - z < 2 Sqrt[y (1 - x - y - z)], {z, y, x}][2]

(where the original constraint is the first argument).

This yielded

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


Using this result as the constraint in a subsequent twofold integration over {x,0,1} followed by {y,0,1} gave an outcome in z, which I expanded out, and then to which I applied term-by-term the function

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


given by Carl Woll in his answer to denestSqrt . The 231 individual terms obtained were then each integrated over z with respect to which of the three ranges ([1/8 (2 - Sqrt2),1/6], [1/6,1/4], [1/4,1/3]) of z within which they fell.

Performing Expand[FullSimplify[Together[]]] on the sum of the resultant 231 integrations yielded (LeafCount=222)

29902415923/497664 - 50274109/(512 Sqrt[2]) - (3619006726645466935245 \[Pi])/(82688 Sqrt[2]) + (926465722021570344469365 ArcCot[Sqrt[2]])/(2646016 Sqrt[2]) - (926465722020908726376075 ArcCot[3 + Sqrt[2]])/(2646016 Sqrt[2]) + (1024176615 ArcTan[2/Sqrt[5]])/(16384 Sqrt[2]) + (48882576431592402695 ArcTan[5/(2 Sqrt[14])])/(194055791689728 Sqrt[2]) - (1861973732913175 Sqrt[2]ArcTan[(4 - Sqrt[2])/Sqrt[18 + 8 Sqrt[2]]])/11844225567 - (48882576431592402695 ArcTan[(20 + 5 Sqrt[2] + 4 Sqrt[7] - 8 Sqrt[14])/Sqrt[1458 + 248 Sqrt[2] + 280 Sqrt[14]]])/(194055791689728 Sqrt[2]) + (1024176615 I Log[1/81 (-46 + 10 Sqrt[10] + I Sqrt[5 (689 + 184 Sqrt[10])])])/(16384 Sqrt[2])


which coincides--using the criterion (suggested by Bob Hanlon in his comment in BlockBlock[{\$MaxExtraPrecision = 1000}, ans1 - ans2 // N[#, 1000] &]--with the much simpler final answer (LeafCount = 36) of JimB.

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

The first two terms of both expressions match.

So, the last eight terms of the former result must be simplifiable to the last two terms of the latter answer.

However, the denestSqrt function--in addition to the Tessore transformation--was not known to me in 2009, so it still remains somewhat of a conundrum as to how that 3D constrained integration was fully carried out. (The result must have also been quite complicated--as indicated by the ensuing correspondence, noted in the statement of the question, with V. Jovovic regarding transformations based on the relation cos(n*phi) =a(n)/3^n, where phi=ArcCos(1/3).)