Compute a certain "separability probability" via a constrained 4D integration over $[-1,1]^4$

Question

I have employed the GenericCylindricalDecomposition command (applied to $6 \times 6$ "density matrix" positivity constraints) to compute the following expression

0 < U < 1 && -1 < a < 1 && -Sqrt[1 - a^2] < b < Sqrt[1 - a^2] && ((-Sqrt[U - b^2 U] < 
  c < -Sqrt[a^2 U] && (a b c)/(-1 + b^2) - 
   Sqrt[(-c^2 + b^2 c^2 + U - 2 b^2 U + b^4 U + a^2 c^2 U - 
    a^2 U^2 + a^2 b^2 U^2)/((-1 + b^2)^2 U)] < 
  d < (a b c)/(-1 + b^2) + 
   Sqrt[(-c^2 + b^2 c^2 + U - 2 b^2 U + b^4 U + a^2 c^2 U - 
    a^2 U^2 + a^2 b^2 U^2)/((-1 + b^2)^2 U)]) || (-Sqrt[a^2 U] < 
  c < Sqrt[
  a^2 U] && (a b c)/(-1 + b^2) - Sqrt[(
   1 - a^2 - 2 b^2 + a^2 b^2 + b^4 - c^2 + a^2 c^2 + 
    b^2 c^2)/(-1 + b^2)^2] < 
  d < (a b c)/(-1 + b^2) + Sqrt[(
   1 - a^2 - 2 b^2 + a^2 b^2 + b^4 - c^2 + a^2 c^2 + 
    b^2 c^2)/(-1 + b^2)^2]) || (Sqrt[a^2 U] < c < Sqrt[
  U - b^2 U] && (a b c)/(-1 + b^2) - 
   Sqrt[(-c^2 + b^2 c^2 + U - 2 b^2 U + b^4 U + a^2 c^2 U - 
    a^2 U^2 + a^2 b^2 U^2)/((-1 + b^2)^2 U)] < 
  d < (a b c)/(-1 + b^2) + 
   Sqrt[(-c^2 + b^2 c^2 + U - 2 b^2 U + b^4 U + a^2 c^2 U - 
    a^2 U^2 + a^2 b^2 U^2)/((-1 + b^2)^2 U)]))

Subject to these constraints, I want to integrate $\frac{3}{2 \pi^2}$ over $a,b,c,d \in [-1,1]$, yielding a function of $U$ (so the Boole command should be employed, it seems). (The $\frac{3}{2 \pi^2}$ comes from enforcement of the positivity of the original $6 \times 6$ density matrix, while the further constraints come from enforcement of the positivity--required for separability--of the "partial transpose" of the density matrix.)

This function when multiplied by

(12 U ((-1 + U) (1 + U (10 + U)) - 6 U (1 + U) Log[U]))/(-1 + U)^5

and integrated over $U \in [0,1]$ should yield a (ten-dimensional) "rebit-retrit Hilbert-Schmidt separability probability" roughly (I guess) equal to 0.6 (cf. sec. IV of https://arxiv.org/abs/1809.09040 ). (An eight-dimensional X-states counterpart problem yields the probability $\frac{16}{3 \pi^2} \approx 0.54038$ https://arxiv.org/abs/1501.02289 p.3.)

If the exact integrations can not be successfully conducted, then it may nevertheless be able to intuit an exact value for the desired probability, by performing a sufficiently high-precision NUMERICAL integration, following the lead of Thies Heidecke in his comment below.

The positivity-based inequalities to which the GenericCylindricalDecomposition command were applied (responding to the request of Heidecke) took the form

1 - a^2 - b^2 > 0 && 1 - a^2 - b^2 > 0 &&
(-1 + a^2) (-1 + c^2) - 2 a b c d - d^2 + b^2 (-1 + d^2) > 0 &&
-c^2 + U - b^2 U > 0 &&
-c^2 - U (-1 + 2 a b c d + d^2 - b^2 (-1 + d^2) + a^2 (-c^2 + U)) > 0 &&
-1 < a < 1 && -1 < b < 1 && -1 < c < 1 && -1 < d < 1 && 0 < U < 1

and after application of FullSimplify

-1 < d < 1 && U < 1 && U > c^2 + b^2 U && a^2 + b^2 < 1 &&
(-1 + a^2) (-1 + c^2) + b^2 (-1 + d^2) > d (2 a b c + d) &&
c^2 + U (-1 + 2 a b c d + d^2 - b^2 (-1 + d^2) + a^2 (-c^2 + U)) < 0  .

Let us omit the constraint

 -c^2 - U (-1 + 2 a b c d + d^2 - b^2 (-1 + d^2) + a^2 (-c^2 + U)) > 0 

in the original (first)--corresponding to the positivity of the determinant of the partial transpose of the $6 \times 6$ density matrix--of these two sets of constraints. (The remaining partial transpose constraint -c^2 + U - b^2 U > 0 corresponds to the positivity of a $5 \times 5$ minor.) Then, based on the associated ("truncated") GenericCylindricalDecomposition,

0 < U < 1 && -1 < b < 1 && -Sqrt[U - b^2 U] < c < Sqrt[U - b^2 U] && -Sqrt[1 - b^2] < a < Sqrt[ 1 - b^2] && (a b c)/(-1 + b^2) - Sqrt[( 1 - a^2 - 2 b^2 + a^2 b^2 + b^4 - c^2 + a^2 c^2 + b^2 c^2)/(-1 + b^2)^2] < d < (a b c)/(-1 + b^2) + Sqrt[( 1 - a^2 - 2 b^2 + a^2 b^2 + b^4 - c^2 + a^2 c^2 + b^2 c^2)/(-1 + b^2)^2]

we can, in fact, proceed analytically, integrating over d, a, c, b, in that order, and obtain an associated "separability function", \begin{equation} \frac{2 \left(\sqrt{(1-U) U}+\sin ^{-1}\left(\sqrt{U}\right)\right)}{\pi }. \end{equation} When multiplied--following our basic algorithm--by

(12 U ((-1 + U) (1 + U (10 + U)) - 6 U (1 + U) Log[U]))/(-1 + U)^5

and integrated over $U \in [0,1]$, we obtain an upper bound of \begin{equation} \frac{919}{5}-264 \log (2) \approx 0.809144 \end{equation} on the separability probability, following the AdaptiveQuasiMonteCarlo approach of Heidecke, that appears to be on the order of 0.744.

Answer 1

Numerically stable integrand

First thing we have to do is make the integrand numerically stable, if we look at the vanilla expression

f = (12 U ((-1 + U) (1 + U (10 + U)) - 6 U (1 + U) Log[U]))/(-1 + U)^5

and its plot

Plot[f, {U, 0, 1}, PlotRange -> {0, 6/5}]

we see that there are no visible singularities, which is good, but $U=0$ might be a problem because of the Log[U] term and $U=1$ might cause problems because of the (-1 + U)^5 part in the denominator. We can fix this by replacing the critical parts with Taylor series approximations

Normal@Series[f, {U, 0, 1}]
FullSimplify@Normal@Series[f, {U, 1, 3}]

12 U

6/35 (5 + U (5 + (-4 + U) U))

which don't contain singularities and combined with the original expression in the numerically uncritical part allow us to construct a smooth version of f

smoothf = Piecewise[
  {{12 U, 0 <= U < 0.00001},
   {f, 0.00001 <= U < 0.97},
   {6/35 (5 + U (5 + (-4 + U) U)), 0.97 <= U <= 1}
  }
] 

Let's check how it compares to the original expression at the critical points

Plot[f - smoothf, {U, 0, 0.000011}, PlotRange -> {-1, 1}/10^6]
Plot[f - smoothf, {U, 0.97, 1}, PlotRange -> {-1, 1}/10^6]

and we see that U==0 probably is uncritical but U==1 we see lots of cancellation error which is not present in our Taylor series and better suitable for integration.

Edit:

Since U==0 is uncritical it's enough to exclude the zero and not do taylor expansion around zero, this will get rid of some error if we want to go to higher precision estimates later. The U==1 part has some drift, too, so we use a higher order exapnsion

FullSimplify@Normal@Series[f, {U, 1, 4}]

1/35 (25 + U (50 + U (-54 + (26 - 5 U) U)))

there to get

smoothf = Piecewise[
  {{f, 0 < U < 0.97},
   {1/35 (25 + U (50 + U (-54 + (26 - 5 U) U))), 0.97 <= U <= 1}
  }
]

Integration strategy

Default

The Boole part of our integrand behaves like a UnitStep with lots of sudden jumps, which are hard to integrate for lots of quadrature schemes, which hope for a somewhat smooth integrand.

We can see that from the warnings when going with the default integration Method

booleregion = Boole[
  U > c^2 + b^2 U && a^2 + b^2 < 1 &&
  (-1 + a^2) (-1 + c^2) + b^2 (-1 + d^2) > d (2 a b c + d) && 
  c^2 + U (-1 + 2 a b c d + d^2 - b^2 (-1 + d^2) + a^2 (-c^2 + U)) < 0
]

NIntegrate[
  3/(2 \[Pi]^2) booleregion smoothf
  , {a, -1, 1}, {b, -1, 1}, {c, -1, 1}, {d, -1, 1}, {U, 0, 1}
  , PrecisionGoal -> 12
] // Timing

NIntegrate::slwcon: Numerical integration converging too slowly; suspect one of the following: singularity, value of the integration is 0, highly oscillatory integrand, or WorkingPrecision too small.

NIntegrate::eincr: The global error of the strategy GlobalAdaptive has increased more than 2000 times. The global error is expected to decrease monotonically after a number of integrand evaluations. Suspect one of the following: the working precision is insufficient for the specified precision goal; the integrand is highly oscillatory or it is not a (piecewise) smooth function; or the true value of the integral is 0. Increasing the value of the GlobalAdaptive option MaxErrorIncreases might lead to a convergent numerical integration. NIntegrate obtained 0.7314353103538273 and 0.07593969354462965 for the integral and error estimates.

and see that it struggled from the discontinuities, but at least we have a first estimate of 0.731 with a quite high error estimate of 0.076.

Monte Carlo

For these kind of non-smooth integrands Monte Carlo integration still guarantees integration error proportional to $n^{-1/2}$ where $n$ is the number of evaluation points regardless of dimension! So let's try that!

Instead of "MonteCarlo" we directly go to "AdaptiveQuasiMonteCarlo" which has at least as good convergence as vanilla Monte Carlo and additionally adaptively subdivides the integration region to focus integrand evaluation mainly on parts where the integrand varies a lot, which makes it even more efficient. With this we get

NIntegrate[
  3/(2 \[Pi]^2) booleregion smoothf
  , {a, -1, 1}, {b, -1, 1}, {c, -1, 1}, {d, -1, 1}, {U, 0, 1}
  , PrecisionGoal -> 12
  , Method -> {"AdaptiveQuasiMonteCarlo", "MaxPoints" -> 10^7}
] // Timing

NIntegrate::maxp: The integral failed to converge after 10000100 integrand evaluations. NIntegrate obtained 0.7440071534377891 and 0.0005174333053340568 for the integral and error estimates.

0.744007

with three digits precision after 2 minutes of computation. Much better! You can increase "MaxPoints" to get higher precision estimates.

Answer 2

We can make progress on this problem by implementing the result--given at the outset of the question--of the GenericCylindricalDecomposition, by integrating over the variables, $d$ and then $c$.

In addition to the three initial inequalities 0 < U < 1 && -1 < a < 1 && -Sqrt[1 - a^2] < b < Sqrt[1 - a^2], the result consists of three parts joined by the OR command ||, each of the three parts consisting of inequalities over $d$ and $c$.

Integrating $\frac{3}{2 \pi^2}$ over $d$, for the first part (as well as the third), we obtain

(3 Sqrt[((-1 + b^2 + a^2 U) (c^2 + (-1 + b^2) U))/((-1 + b^2)^2 U)])/\[Pi]^2,

while for the second, we find

(3 Sqrt[((-1 + a^2 + b^2) (-1 + b^2 + c^2))/(-1 + b^2)^2])/\[Pi]^2 .

Now, integrating--subject to the first three inequalities of the GenericCylindricalDecomposition--these outcomes over $c$, we obtain for the first of the two integrations

$A=$ \begin{equation} \frac{3 \left( \begin{array}{cc} \{ & \begin{array}{cc} \frac{\sqrt{-U \left(U a^2+b^2-1\right)} \left(-2 \tan ^{-1}\left(\frac{a}{\sqrt{-a^2-b^2+1}}\right) \left(b^2-1\right)+\pi \left(b^2-1\right)+2 a \sqrt{-a^2-b^2+1}\right)}{2 \left(b^2-1\right)} & a\geq 0 \\ \frac{\sqrt{-U \left(U a^2+b^2-1\right)} \left(2 \tan ^{-1}\left(\frac{a}{\sqrt{-a^2-b^2+1}}\right) \left(b^2-1\right)+\pi \left(b^2-1\right)-2 a \sqrt{-a^2-b^2+1}\right)}{2 \left(b^2-1\right)} & \text{True} \\ \end{array} \\ \end{array} \right)}{2 \pi ^2} , \end{equation}

while for the second, we find

B = \begin{equation} \text{ConditionalExpression}\left[\frac{3 \sqrt{-a^2-b^2+1} \left(-2 \left| a\right| \sqrt{-U \left(a^2 U+b^2-1\right)}+\left(b^2-1\right) \left(2 i \log \left(\sqrt{U} \left| a\right| +i \sqrt{a^2 (-U)-b^2+1}\right)-i \log \left(1-b^2\right)+\pi \right)\right)}{2 \pi ^2 \left(b^2-1\right)},a\neq 0\right]. \end{equation}

To derive the desired separability probability, we must now integrate over $U,a,b$ (employing again the first three inequalities of the GenericCylindricalDecomposition), the product of $2 A$, and of $B$, with

C= (12 U ((-1 + U) (1 + U (10 + U)) - 6 U (1 + U) Log[U]))/(-1 + U)^5 .

(We bypass for the moment the search for the "separability function" that would be obtained by further integrating $2 A$ and $B$--without taking their products with $C$--over $b$ and $a$, yielding a function of $U$.)

We have so far been unable to perform these integrations exactly. Numerical integrations clearly indicate (quite interestingly) that the result of integrating $2 A C$ is the same as of integrating $B C$--so, we can use either $4 A C$ or $2 B C$ in the total computations. A WorkingPrecision->20 (4-minute) calculation gave us for the separability probability the result 0.74792453608226129700,

NIntegrate::eincr: The global error of the strategy GlobalAdaptive has increased more than 2000 times. The global error is expected to decrease monotonically after a number of integrand evaluations. Suspect one of the following: the working precision is insufficient for the specified precision goal; the integrand is highly oscillatory or it is not a (piecewise) smooth function; or the true value of the integral is 0. Increasing the value of the GlobalAdaptive option MaxErrorIncreases might lead to a convergent numerical integration. NIntegrate obtained 0.7479245360822612969975167328445174396108540562661672111470100755943716 and 2.751458747155140755116865153269220994697828764362273923451154625009759*10^-10 for the integral and error estimates.

while the (2-minute) three-digit-precision result of Heidecke in his answer to the question was 0.744007.

We still aspire to exact integration results, and if not so, the possible discernment of an underlying exact formula for the separability probability--in correspondence with the $\frac{919}{5}−264\log(2)≈0.809144$ upper bound given near the end of the original question, and the associated separability function \begin{equation} s(U)= \frac{2 \left(\sqrt{(1-U) U}+\sin ^{-1}\left(\sqrt{U}\right)\right)}{\pi }. \end{equation}

Let us note that for $U=1$, this function yields $\sqrt{\frac{2}{\pi}} \approx 0.797885$, while numerics strongly indicate that for the separability probability function ($S(U)$)--not yet constructed--based on the full set of partial transpose positivity constraints, the corresponding $U=1$ value is simply $\frac{1}{2}$. Also, $\int_{U=0}^1 s(1) dU= \frac{3}{4}$.

Further examination and analysis--to be expanded upon--leads us to the conclusion that the formula for the titular separability probability is \begin{equation} -\frac{34 C}{3}+\frac{5609}{30}-264 \log (2)+\frac{1}{6} \pi (7-9 \pi -97 \log (2)+93 \log (3)) \approx 0.747924541771417135649643, \end{equation} where $C$ is Catalan's constant $\approx 0.915966$.

Answer 3

Thanks to the skillful answer of "Student" to the question posed in https://mathoverflow.net/questions/325697/compute-via-possibly-integration-by-parts-the-definite-integral-of-the-product#325715

we are able to, first, conclude that the desired "separability function"--incorporating all the positivity constraints (on the $6 \times 6$ density matrix and its "partial transpose")--is

A = -((2 ((-1 + U) (-Sqrt[U] + (1 + U) ArcTanh[Sqrt[U]]) - 4 U PolyLog[2, Sqrt[U]] + U PolyLog[2, U]))/([Pi]^2 U)).

(Let us note that with the transformation $U \rightarrow \varepsilon^2$, this becomes identically equal to the "two-rebit" function $\tilde{\chi}_1(\varepsilon)$ of Lovas and Andai, reported as eq. (2) in https://arxiv.org/abs/1701.01973 .)

Then, the integration over $U \in [0,1]$ of the product of $A$ with the previously noted

B= (12 U ((-1 + U) (1 + U (10 + U)) - 6 U (1 + U) Log[U]))/(-1 + U)^5

gives us the desired "two-rebit separability probability".

Student was able to demonstate (working with a $U \rightarrow u^2$ transformation) that this probability was equal to \begin{equation} \frac{256-168 \zeta (3)}{\pi ^2}-38+48 \log (2) \approx 0.747925. \end{equation}

Further, seeing the relevance here of the $\tilde{\chi}_1(\varepsilon)$ function of Lovas and Andai, it seems not too far a leap for us to advance hypotheses that--based on the two-qubit function $\tilde{\chi}_2(\varepsilon)=\frac{1}{3} \varepsilon^2 (4 -\varepsilon^2)$, the "two-qubit separability probability" is $\frac{5}{3} \left(112 \pi ^2-1105\right) = \frac{560 \pi ^2}{3}-\frac{5525}{3} \approx 0.659488$, and, further, with $\tilde{\chi}_{4}(\varepsilon)=\frac{1}{35} \varepsilon ^4 \left(15 \varepsilon ^4-64 \varepsilon ^2+84\right)$, the "two-quater[nionic]-bit separability probability" is $\frac{8962661573}{4725}-192192 \pi ^2 \approx 0.583115$. (We observe that $1105 = 5 \cdot 13 \cdot 17$, and $112=2^4 \cdot 7$, while $192192-2^6 \cdot 3 \cdot 7 \cdot 11 \cdot 13$, and $4725 = 3^3 \cdot 5^2 \cdot 7$.)

We are presently investigating this two-qubit hypothesis, but it presently seems quite challenging both in numerical and analytic terms.