Integration of a function over 7-sphere Suppose we have  $x_1^2 + y_1^2 + x_2^2 + y_2^2 + x_3^2 + y_3^2 + x_4^2 + y_4^2 = 1$  and we define  $z_j = x_j + iy_j$,  where  $j = 1,\,2,\,3,\,4$.
The problem is finding or approximating the following integral (The actual problem is more complex than this !)
$$I(k) = \int\limits_{\mathbb{S}^7}\lvert(z_1z_4 - z_2z_3)^k\rvert,$$ where $\mathbb{S}^7$ is the $7$-sphere, and the integral is over its surface measure. The ultimate aim is to find or get an estimate of the ratios of the type  $\displaystyle \frac{I(k)}{I(k+1)}$.
I have never done an integration over such a surface, the hypersphere and don't know where to begin. Any hint or any pointers would really help me to get started.
(I have been suggested to make use of the fact that the points $x_i$'s and $y_i$'s lie on a sphere; hence we can use it to reduce one variable and make use of Fubini's Theorem with proper limits for integration to proceed. But I can't see how to make use of this.)

What have I done till now?
As the problem say, it would be acceptable if we can approximate the ratios or the integral, so
$$\lvert z_1z_4 - z_2z_3\rvert \, \leq \, \lvert z_1z_4\rvert + \lvert z_2z_3\rvert \, = \, \lvert z_1\rvert\lvert z_4\rvert + \lvert z_2\lvert\rvert z_3\rvert.$$
This converts everything into real variables, and then I have used Mathematica to calculate this. But the problem is Mathematica gives the solution for small values of  $k$, but as I have to find it for the general case, it is not of much use.
Also, I have tried using the spherical coordinates to help with the integral, but they make it more complicated.
Any pointers on this front are also welcome.

Here are the values of $I(k)$  for  $k=1,\,\cdots,\,10$.
$z_1 = x_1 + ix_2, \ \ z_2 = x_3+ix_4, \ \ z_3 = x_5+ix_6, \ \ z_4 = x_7+ix_8$  then
$$\mid z_1z_4 - z_2z_3| \ \ = \ \ \sqrt{\left(x_4 x_5+x_3 x_6-x_2 x_7-x_1 x_8\right){}^2+\left(x_3 x_5-x_4 x_6-x_1 x_7+x_2 x_8\right)^2} $$
To integrate over $\mathbb{S}^7$, I have used two methods, one is using NIntegrate and the second using integrateSphere function from HFT Software package.
(Note : HFT package doesn't work for odd values of $k$, but for even values of $k$, the results match perfectly with that of INtegrate)
$$\begin{align}
I(1) \ \ &\approx \ \ 9.54903\\
I(2) \ \ &= \ \ \frac{\pi^4}{3\times 10} \ \ \approx \ \ 3.24697\\\\
I(3) \ \ &\approx \ \ 1.19691\\
I(4) \ \ &= \ \ \frac{\pi^4}{3\times70} \ \ \approx \ \ 0.463853\\\\
I(5) \ \ &\approx \ \ 0.188085\\
I(6) \ \ &= \ \ \frac{\pi^4}{3\times420} \ \ \approx \ \ 0.0773088\\\\
I(7) \ \ &\approx \ \ 0.0329181\\
I(8) \ \ &= \ \ \frac{\pi^4}{3\times2310} \ \ \approx \ \ 0.0140561\\\\
I(9) \ \ &\approx \ \ 0.00620447\\
I(10) \ \ &= \ \ \frac{\pi^4}{3\times12012} \ \ \approx \ \ 0.0027031\\\\
\end{align}$$
 A: $\newcommand{\R}{\mathbb{R}}\newcommand{\C}{\mathbb{C}}$Here's an approach I like that takes advantage of the fact that the integrands are homogeneous. I'll just show how to compute
$$ \int_{S^7} |z_1z_4|. $$
First, denote $\hat{z} \in \C^4$ and $|\hat z|^2 = |\hat z_1|^2 + |\hat z_2|^2 + |\hat z_3|^2 + |\hat z_4|^2$, then, using polar coordinates
\begin{align*}
\int_{\R^8} e^{-|\hat z|^2}|\hat z_1\hat z_4|\,d\hat{z} 
&= \int_{z \in S^7}\int_{r=0}^{r=\infty} e^{-r^2}r^2|z_1z_4|r^{7}\,dr\,dz\\
&= \left(\int_{r=0}^{r=\infty} r^9e^{-r^2}\,dr\right)\int_{S^7} |z_1z_4|\,dz,
\end{align*}
where $dz$ is the surface area measure on $S^7$. The second factor on the right is what we want to calculate. The value of the first factor on the right is easy to evaluate in terms of the gamma function. So it remains only to evaluate the left integral. But this is straightforward. For convenience, I'll write $z \in \C^4$ instead of $\hat{z}$. Then
\begin{align*}
\int e^{-|z|^2}|z_1z_4|\,dz
&= \int_{\C} e^{-|z_1|^2}|z_1|\,dz_1\int_{\C} e^{-|z_4|^2}|z_4|\,dz_4 \int_{\C} e^{-|z_2|^2}\,dz_2\int_\C e^{-|z_3|^2}\,dz_3\\
&= \left(\int_{\C} e^{-|z|^2}|z|\,dz\right)^2\left(\int_{\C} e^{-|z|^2}\,dz\right)^2.
\end{align*}
The two integrals on the right are also easily evaluated.
A: This is not a complete answer but shows how to simplify your integral:
As pointed out by mlk in the comments the Hopf fibration plays a crucial role: Consider $\mathbb S^7\subset\mathbb C^4.$ On $\mathbb C^4=\mathbb C^2\oplus\mathbb C^2$ the special unitary group $\mathrm{SU}(2)$ acts by the direct sum of the standard representation. Note that $\mathrm{SU}(2)=\mathbb S^3.$ The Hopf fibration $$\mathrm{SU}(2)\to\mathbb S^7\to\mathbb S^4$$ is given by taking the quotient of this action restricted to $\mathbb S^7\subset\mathbb C^4.$ Then, one directly computes that $$z_1z_4-z_2z_3$$ is invariant under the action of $\mathrm{SU}(2).$ Hence, by Fubini using the fact that we consider a Riemannian submersion, the integral $I(k)$ is given by a constant factor $C\neq0$ (comming from the volume of $\mathrm{SU}(2)$ ) times an integral over the quotient $\mathbb S^4=\mathbb S^7/\mathrm{SU}(2).$  For your final question concernig $I(k)/I(k+1)$, the constant $C$ does not play any role.
Let us identify the 4-sphere with the quaterninic projective line $\mathbb H P^1.$ Using affine coordinates (being basically equivalent to  coordinates using a suitable stereographic projection, your integrand becomes $$|z_1 z_2-z_2 z_3|^k( [x_1+x_2i+(x_3+x_4 i) j,1])=|x_3+x_4 i|^k.$$ The volume form (what you call the surface measure) on the quotient is the standard volume form and can be written as $$vol= \frac{4}{(1+x_1^2+x_2^2+x_3^2+x_4^2)^4} dx_1\wedge dx_2\wedge dx_3\wedge dx_4.$$ So you should integrate
$$\int_{\mathbb R^4}\frac{4 |x_3+x_4 i|^k}{(1+x_1^2+x_2^2+x_3^2+x_4^2)^4} dx_1\wedge dx_2\wedge dx_3\wedge dx_4.$$ I do not have time to finish the computation, but I am sure the last integral can be simplified further.
A: The result is
$$
I(k)=\frac{2\pi^4\ \Gamma\left(\frac{k}{2}+1\right)\Gamma\left(\frac{k}{2}+2\right)}{\Gamma(k+4)}\ ,
$$
which can be simplified a bit more using the Legendre Duplication Formula.
More generally, we have
$$
J(n,k):=\int_{M_n(\mathbb{C})}e^{-\sum_{i,j}|z_{ij}|^2}|{\rm det}((z_{ij})_{1\le i,j\le n})|^{k}
\prod_{i,j}\frac{d\Re z_{ij}d\Im z_{ij}}{\pi}
=\prod_{i=1}^{n}\frac{\Gamma\left(\frac{k}{2}+i\right)}{\Gamma(i)}\ .
$$
Specializing to $n=2$ and going to (hyper)spherical coordinates, we have
$$
J(2,k)=\frac{1}{\pi^4}I(k)\times\int\limits_{0}^{\infty} e^{-r^2}r^{2k+7}dr=\frac{1}{2\pi^4}I(k)\Gamma(k+4)\ ,
$$
from which the formula for $I(k)$ easily follows.
Remark:
A very useful reference for explicit computations similar to this MO question (e.g., emphasizing the use of Carlson's Theorem to reduce to the case where $k$ is an even integer, and then Isserlis-Wick Feynman diagrammatic computations) is the article "Hauteur des hypersurfaces et fonctions Zêta d'Igusa" by Julien Cassaigne and Vincent Maillot. I first learned the Carlson$\rightarrow$combinatorics trick from their paper. However, I just saw in the article "The importance of the Selberg integral" by Forrester and Warnaar that Selberg used the same idea.
Finally, note that the 2nd displayed formula on page 95 of the book "An Introduction to the Theory of Local Zeta Functions" by Jun-Ichi Igusa for the quantity $J(n,k)$ is not correct.
Having a formula for $J(n,k)$ is related to the explicit knowledge of the relevant Bernstein-Sato polynomial, because of the so-called Cayley identity
$$
{\rm det}(\partial Z)\ [{\rm det}(Z)]^k=k(k+1)\cdots(k+n-1)\ [{\rm det}(Z)]^{k-1}\ .
$$
For a wealth of information on this identity see the article "Algebraic/combinatorial proofs of Cayley-type identities for derivatives of determinants and pfaffians" by Caracciolo, Sokal and Sportiello.

Addendum:
To make this more self-contained, here is a very general method to compute $J(n,k)$ (and a whole lot more).
Let $z^{(1)}\ldots,z^{(n)}$ denote the row vectors of the matrix $Z=(z_{ij})_{i,j}$.
Suppose $k=2p$ with $p$ a nonnegative integer.
We have to compute a Gaussian expectation of the integrand $|{\rm det}(Z)|^{2p}=({\rm det}(ZZ^{\ast}))^p$. Note that
$$
{\rm det}(ZZ^{\ast})=\sum_{\sigma\in S_n}{\rm sgn}(\sigma)
\langle z^{(1)},z^{(\sigma(1))}\rangle\cdots\langle z^{(n)},z^{(\sigma(n))}\rangle
$$
and thus involves (unnormalized) antisymmetrization. The integrand is the previous expression to the $p$-th power, so we have $p$ antisymmetrizers.
The Isserlis-Wick (IW) formula for complex Gaussian measures can be found for instance in
my article "Feynman Diagrams in Algebraic Combinatorics" (Theorem 1).
For each of the independent random Gaussian vectors $z^{(i)}$, we have as a consequence of the IW formula
$$
\mathbb{E}\ z^{(i)}_{j_1}\cdots z^{(i)}_{j_p}\bar{z}^{(i)}_{\ell_1}\cdots \bar{z}^{(i)}_{\ell_p} 
=\sum_{\tau\in S_p} \delta_{j_1,\ell_{\tau(1)}}\cdots \delta_{j_p,\ell_{\tau(p)}}\ ,
$$
where I used Kronecker deltas.
So now we also have $n$ (unnormalized) symmetrizers.
For the last step, I will need the reader to be an aficionado of diagrammatic algebra. Using pictures I will not draw here, it is trivial to see that what we are computing is the trace of the linear map $(\mathbb{C}^{n})^{\otimes np}\rightarrow (\mathbb{C}^{n})^{\otimes np}$ corresponding (by the so-called Schur-Weyl duality) to the group algebra element in $\mathbb{C}[S_{np}]$ given by a Young symmetrizer $Y(T)$ of a Young tableau $T$ of rectangular shape $\lambda=(p^n)$.
This trace is given by
$$
\prod_{\square\in T}(n+c(\square))\ ,
$$
where $c(\square)$ is the content of the box $\square$, namely, $j-i$ if $\square$ is in the $i$-th row and $j$-th column of the diagram.
This follows from the hook-content formula. A very nice proof is in Appendix B.4 of the book "Group Theory" by Cvitanović.
We now have
$$
J(n,k)=\prod_{\square\in T}(n+c(\square))\ .
$$
Collecting the factors row-wise from top to bottom, for the shape $\lambda=(p^n)$, we get
$$
J(n,k)=\frac{(n+p-1)!}{(n-1)!}\times\frac{(n+p-2)!}{(n-2)!}\times\cdots\times\frac{p!}{0!}\ ,
$$
i.e.,
$$
J(n,k)=\prod_{i=1}^{n}\frac{\Gamma\left(\frac{k}{2}+i\right)}{\Gamma(i)}\ .
$$
The above only uses the rectangular partition case of the hook-content formula, but it should work for any shape. The same method for general partitions should work for the following integral
$$
\int_{M_{n,N}(\mathbb{C})}e^{-{\rm tr}(ZZ^{\ast})}
\Delta_1^{\lambda_1-\lambda_2}\Delta_2^{\lambda_2-\lambda_3}\cdots
\Delta_{n-1}^{\lambda_{n-1}-\lambda_n}\Delta_n^{\lambda_n}
\ \prod_{i=1}^{n}\prod_{j=1}^{N}\frac{d\Re z_{ij}d\Im z_{ij}}{\pi}
$$
over matrices $Z$ of format $n\times N$ with $N\ge n$. Here $\Delta_k$ denotes the minor determinant of size $k$ in the top left corner of the matrix $ZZ^{\ast}$. Namely, the result should be
$$
\prod_{j=1}^{n}\frac{\Gamma(N+\lambda_{j}-j+1)}{\Gamma(N-j+1)}\ .
$$
Since the crux of the method is the use of the hook-content formula, i.e., the evaluation of $s_{\lambda}(1,\ldots,1)$, I suspect this can be generalized to a $q$-deformed version thanks to Stanley's hook-content formula, i.e., the evaluation of $s_{\lambda}(1,q,\ldots,q^{N-1})$.
