Speeding up trigonometric integral

Question

Context

On a possible non trivial toric topology for the Universe (nothing less!).

Problem

I would like to carry out the following integral for $\ell=2,4\cdots 20$.

$$\int _0^{\pi }\int _0^{2 \pi } \sin(\theta ) P_{\ell }\left(\frac{1}{2} \sin(\theta ) (\cos(\phi )-\sin(\phi ))\right) P_{\ell }\left(\frac{1}{2} \sin(\theta ) (\cos(\phi )+\sin (\phi ))\right) \times $$ $$P_{\ell }\left(\frac{\sin (\theta ) ((K+1) \cos (\phi )+(1-K) \sin (\phi ))}{2 \sqrt{K^2+2 \mu K+1}}\right)d\phi d\theta$$ where $P_\ell$ are Legendre Polynomials.

For $\ell$ larger than say 8, Mathematica takes forever and runs out of memory. I have found out a method to circumvent the memory problem (see attempt below) but it still takes very long time to carry out the integration for $\ell>10$.

Question

Is there a way to be smart about this class of integral? Another approach to the one below?

Attempt

I have defined the integrant as

integ[ℓ_] := 
 LegendreP[ℓ, (Cos[p] - Sin[p]) Sin[t]/
     2] LegendreP[ℓ, (Cos[p] + Sin[p]) Sin[t]/
     2] LegendreP[ℓ, ((1 + K) Cos[p] + (1 - K) Sin[p]) Sin[
       t]/2/Sqrt[1 + 2 K μ + K^2]] Sin[t]

And the following integration rule:

r1 = {Exp[ Complex[0, b_] t +  p Complex[0, c_] ] -> 
   Integrate[Exp[ I b  t + I c p], {t, 0, Pi}, {a, 0, 2 Pi}],
  Exp[ Complex[0, b_] t ] -> 
   Integrate[Exp[ I b  t], {t, 0, Pi}, {a, 0, 2 Pi}],
  Exp[   p Complex[0, c_] ] -> 
   Integrate[Exp[ I c p], {t, 0, Pi}, {a, 0, 2 Pi}]}

so that I expand my integrant into sines and cosines into complex exponentials, and carry out the integration via substitution:

   integ[2] // TrigToExp // Expand // 
       Collect[#, Exp[_]] & // (# /. r1) & // Apart // Simplify

 (*-((Pi*(43*K^2 + 84*K*μ + 43))/
  (560*(K^2 + 2*K*μ + 1))) *)

A simple timing

Table[{i, integ[i] // TrigToExp // Expand // 
                 Collect[#, Exp[_]] & // (# /. r1) & // Apart // 
   Simplify;//Timing}, {i, 2, 6,2}]

suggest a $n^{4.5}$ scaling.

Answer 1

Use the following representation of the Legendre polynomials: $$ P_n(x) = 2^n \sum_{k=0}^n x^k \binom{n}{k} \binom{\frac{n+k-1}{n}}{n} $$ Note that the sum effectively is over $k \equiv n \bmod 2$.

Expand each Legendre polynomial into a sum. Integration with respect to $\theta$ is easy: $$ \int_0^{\pi} \sin^{k_1+k_2+k_3+1} \theta \mathrm{d}\theta = \operatorname{Beta}\left(\frac{k_1+k_2+k_3+2}{2}, \frac{1}{2}\right) \tag{1} $$ Integration with respect to $\phi$ is more involved. We need the following result, for $p,q \in \mathbb{Z}_{\geqslant 0}$: $$ \int_0^{2\pi} \cos^{p} \phi \sin^q \phi \mathrm{d}\phi = 2 \cos^2\left( \frac{\pi p}{2} \right) \cos^2\left( \frac{\pi q}{2} \right) \operatorname{Beta}\left(\frac{1+p}{2}, \frac{1+q}{2} \right) \tag{2} $$ notice that the right hand side vanishes whenever $p$ or $q$ are odd.

We now apply binomial theorem to the following and use eq. $(2)$ $$ \int_0^{2\pi} \left(\cos \phi - \sin\phi\right)^{k_1} \left(\cos \phi + \sin\phi\right)^{k_2} \left((1+\kappa)\cos \phi + (1-\kappa)\sin\phi\right)^{k_3} \mathrm{d}\phi $$ Combining these results in the following Sum representation on your integral:

binoms[n_, k_] := Binomial[n, k] Binomial[(n + k - 1)/2, n];

sumint[n_Integer, kappa_, mu_] := 
 Sum[If[EvenQ[k1 + k2 + k3], 
   binoms[n, k1] binoms[n, k2] binoms[n, 
     k3] Beta[(k1 + k2 + k3)/2 + 1, 1/2] 2^(3 n + 1 - k1 - k2 - k3)
     1/(1 + 2 kappa mu + kappa^2)^(k3/2)
     If[EvenQ[m1 + m2 + m3], 
     Binomial[k1, m1] Binomial[k2, m2] Binomial[k3, m3] (-1)^
      m1 (1 + kappa)^(k3 - m3) (1 - kappa)^
      m3 Beta[(1 + k1 - m1 + k2 - m2 + k3 - m3)/2, (1 + m1 + m2 + m3)/
       2], 0], 0], {k1, Mod[n, 2], n, 2}, {k2, Mod[n, 2], n, 2}, {k3, 
   Mod[n, 1], n, 2}, {m1, 0, k1}, {m2, 0, k2}, {m3, 0, k3}, 
  Method -> "Procedural"]

Verification:

integ[n_, kappa_, mu_] := 
 LegendreP[n, 1/2 (Cos[p] - Sin[p]) Sin[t]] LegendreP[n, 
   1/2 (Cos[p] + Sin[p]) Sin[t]] LegendreP[
   n, (((1 + kappa) Cos[p] + (1 - kappa) Sin[p]) Sin[t])/(
   2 Sqrt[1 + 2 kappa mu + kappa^2])] Sin[t]

First check against Integrate for low values of n:

In[8]:= mapint[e_Plus] := 
  Integrate[#, {t, 0, Pi}, {p, 0, 2 Pi}] & /@ e;
mapint[e_] := Integrate[e, {t, 0, Pi}, {p, 0, 2 Pi}]

In[10]:= Table[
   sumint[n, ka, mu] - 
    mapint[Expand[integ[n, ka, mu], _Sin | _Cos]], {n, 0, 3}] // 
  Simplify // AbsoluteTiming

Out[10]= {25.681514, {0, 0, 0, 0}}

Here are timings:

In[45]:= Table[{n, First[AbsoluteTiming[sumint[n, ka, mu];]]}, {n, 1, 
  12}]

Out[45]= {{1, 0.}, {2, 0.}, {3, 0.010000}, {4, 0.050001}, {5, 
  0.070001}, {6, 0.270005}, {7, 0.420008}, {8, 1.180024}, {9, 
  1.510030}, {10, 3.510070}, {11, 4.340087}, {12, 9.130183}}

Evaluation for $n=20$ takes 2 minutes on a real fast machine:

In[14]:= AbsoluteTiming[r20 = sumint[20, ka, mu]; ]

Out[14]= {120.752106, Null}

In[15]:= AbsoluteTiming[r20s=Simplify[r20];]

Out[15]= {0.148191, Null}

In[16]:= LeafCount[r20s]

Out[16]= 385

The expression is not fit for display here, as coefficients are some 20-digit integers. I can post the result if requested.