NIntegrate with a complicated function over the 4-sphere

Question

I would like to maximize the following expression: $$\sum_{k = 1}^4\int_{\mathbb{S^4}}|u - \mathbf{e}_k.u_k||\nabla u_k|.\Big(1 - \frac{1}{8}\Big(\langle{v, \mathbf{e}_k}\rangle^2 + \langle{ e_k,\frac{\nabla u_k}{|\nabla u_k|}}\rangle^2 + \langle{\frac{\nabla u_k}{|\nabla u_k|}, \frac{u - u_k.\mathbf{e}_k}{u - u_k.\mathbf{e}_k}}\rangle^2 + \langle{\frac{u - u_k.\mathbf{e}_k}{u - u_k.\mathbf{e}_k}, v}\rangle^2\Big)\Big).dS$$ where $dS$ denotes the normalised measure of $\mathbb{S}^4$, $u = u_1.\mathbf{e}_1 + u_2.\mathbf{e}_2 + u_3.\mathbf{e}_3 + u_4.\mathbf{e}_4$, where $\mathbf{e}_i$ are the standard vector and $\nabla$ is the gradient on $\mathbb{S}^4$. For the maximization we ask that $|u|^2 = u_1^2 + u_2^2 + u_3^2 + u_4^2$ integrates to $1$ over $\mathbb{S}^4$ AND that $u_i$ are spherical harmonics of degree 3 over $\mathbb{S}^4$.

To do this with Mathematica, I am using a package by Sheldon Axler (https://www.axler.net/HFT_Math.html) that generates an orthonormal basis of spherical harmonics of degree 3 for me. The first step is to integrate the term for $k = 1$.

m = 4;
n = 3;
setDimension[x, n];
d = dimHarmonic[m, n];
basis = basisH[m, x, 
   Sphere] /. \
{\[LeftDoubleBracketingBar]x\[RightDoubleBracketingBar]^2 -> 
    x^2 + y^2 + z^2 + t^2, x[1] -> x, x[2] -> y, x[3] -> z, x[4] -> t};
weight[a_?(VectorQ[#, NumericQ] &), b_?(VectorQ[#, NumericQ] &), 
   c_?(VectorQ[#, NumericQ] &), d_?(VectorQ[#, NumericQ] &)] := 
  1 - 1/8*(x^2 
      + ({x, y, z, t} . unit[{0, b . basis, c . basis, d . basis}])^2
      + (unit[spheregrad[a . basis]] . {1, 0, 0, 0})^2
      + (unit[spheregrad[a . basis]] . 
         unit[{0, b . basis, c . basis, d . basis}])^2);
term[a_?(VectorQ[#, NumericQ] &), b_?(VectorQ[#, NumericQ] &), 
   c_?(VectorQ[#, NumericQ] &), 
   d_?(VectorQ[#, NumericQ] &)] := (1/surfaceArea[n])*
   NIntegrate[
    Sqrt[(b . basis)^2 + (c . basis)^2 + (d . basis)^2]*
     norm[spheregrad[a . basis]]*
     weighta[a, b, c, d], {x, y, z, t} \[Element] Sphere[n], 
    WorkingPrecision -> 15, MaxPoints -> 10000];

Here unit denotes the re-scaling of a vector to norm $1$, spheregrad is the function that return the spherical gradient. Evaluating the function term at a random point on $\mathbb{S}^{64}$ ($d = 16$ is the dimension of spherical harmonics of degree $3$ on $\mathbb{S}^4$, so writing $(a,b,c,d) \in \mathbb{S}^{64}$ with $u_1 = a.basis$, $u_2 = b.basis$ and so on) simply gives Warnings `` and the function returns no answer. I suspect the reason is that the function weight is complicated. Any help in making term return some value is much appreciated!!

Answer 1

As mentioned in the comments, hyperspherical coordinates give a direct parametrization of the surface of the unit 4-sphere. To integrate a function $f(x,y,z,w)$ on this hyperspherical surface $x^2+y^2+z^2+w^2=1$, we use the coordinates

$$ x = \cos (\theta )\\ y = \sin (\theta ) \cos (\chi )\\ z = \sin (\theta ) \sin (\chi ) \cos (\phi)\\ w = \sin (\theta ) \sin (\chi ) \sin (\phi ) $$

to give the 4-surface integral

NIntegrate[f[Cos[θ], Sin[θ]*Cos[χ], Sin[θ]*Sin[χ]*Cos[φ], Sin[θ]*Sin[χ]*Sin[φ]] *
           Sin[θ]^2*Sin[χ], {θ, 0, π}, {χ, 0, π}, {φ, 0, 2π}]

To get more uniform sampling, use $(x,y,\phi)$ as parameters and thus eliminate the Jacobian:

NIntegrate[f[x, y, Sqrt[1-x^2-y^2]*Cos[φ], Sqrt[1-x^2-y^2]*Sin[φ]],
           {x, -1, 1}, {y, -Sqrt[1-x^2], Sqrt[1-x^2]}, {φ, 0, 2π}]

There is a lot of symmetry in this integral, and you are allowed to permute all arguments if it makes your integral simpler.