Taking the surface integral $\int_{S_1(0)}y_jy_k \ d\sigma(y)$ I'm learning how to take surface integrals on the surface of spheres in $\mathbb{R}^n$. This question is related to Calculating the surface integral $\int_{S_1(0)}y_j \ d\sigma(y)$ where I try to compute a similar integral but I don't know if I did everything ok. (read that answer to have a better understanding of how to integrate over a surface if you need)
I'm asked to compute
$$\int_{S_1(0)}y_jy_k \ d\sigma(y)$$ 
Let's imagine the north hemisphere first, through the parametrization $\Sigma(y) = \left(y,\sqrt{1-|y|^2}\right)$:
$$\int_{S_1(0)^+}y_jy_k \ d\sigma(y) = \int_{C_1(0)+}y_jy_k\det [ye_1\ \cdots \  ye_{n-1} \ n(y)]\ dy$$ 
where $C_1(0)$ is just the 'circunference' on which the 'sphere' is parametrized. We can think of it as the region $|y|<1$ for $y\in\mathbb{R}^{n-1}$
Remember that the unit normal $n(y)$ in the unit sphere is just $y$. So the determinant above would give us $y_1\cdots y_{n-1}\sqrt{1-(y_1+\cdots y_{n-1})^2}$ because these are the diagonal terms multiplied and the rest of the elements are $0$ (except for the normal column but it gets multiplied by the other $0$s)
So we end up with
$$\int_{S_1(0)^+}y_jy_k \ d\sigma(y) = \int_{C_1(0)+} y_1\cdots y_j^2 \cdots y_k^2 \cdots \ y_{n-1}\sqrt{1-(y_1^2+\cdots +  y_{n-1}^2)}\  dy$$
Now, breaking these onto the region $C_1(0)$ we get:
$$\int_{S_1(0)^+}y_jy_k \ d\sigma(y) = \\ \int_{-1}^1\cdots\int_{-1}^1 y_1\cdots y_j^2\cdots y_k^2 \cdots \ y_{n-1}\sqrt{1-(y_1^2+\cdots +y_{n-1}^2)} \ dy_1\cdots dy_j \cdots dy_{n-1}$$
And for the south hemisphere:
$$\int_{S_1(0)^-}y_jy_k \ d\sigma(y) =\\ \int_{-1}^1\cdots\int_{-1}^1 y_1\cdots y_j^2 \cdots y_k^2 \cdots \ y_{n-1}\sqrt{-1+(y_1^2+\cdots + y_{n-1}^2)} \ dy_1\cdots dy_j \cdots dy_{n-1}$$
I think it helps if I do the two cases of integration for the north hemisphere first. If $k\neq j$ then I have to integrate $y_i\sqrt(\cdots)$ or $y_i^2\sqrt(\cdots)$. If $k=j$ then it's a matter of integrating $y_i\sqrt(\cdots)$ and $y_i^3\sqrt(\cdots)$. 
So for $k\neq j$, if we integrate $y_i\sqrt(\cdots)$:
$$\frac{1}{2}\int_{-1}^12y_i\sqrt{1-(y_1^2 + \cdots + y_i^2 + \cdots + y_{n-1}^2)}\ dy_i = \\\frac{1}{2}\int_{1}^{1}\sqrt{1-(y_1^2 + \cdots + u + \cdots + y_{n-1}^2)}\ du$$
I'm integrating with the substitution $u = y_i^2$ from $u(-1) = 1$ to $u(1) = 1$ so it should be $0$. This gets multiplied with every other integral so the entire integral is $0$ which is ok according to my book.
Now for $k=j$, if we integrate $y_i^4\sqrt(\cdots)$ we would get something, but it would also be multiplied by the integral of $y_i\sqrt(\cdots)$, so it should be $0$ too. But my book says it is $n^{-1}\int{S_1(0)}1$.
Batominovski's Comment: In the paragraph above, the claim that "it would also be multiplied by the integral of $y_i\sqrt(\cdots)$" is wrong.  We do not have $\int\,(fg) =\left(\int\,f\right)\,\left(\int\,g\right)$.  So, $\int\,g=0$ does not imply $\int\,(fg)=0$.
What is wrong?
 A: You have made some errors in setting up the integral, both here and in your previous post. The parameterization of the upper hemisphere we're using is
$$\Sigma (y) = (y,(1-|y|^2)^{1/2},$$
for $y$ in the open unit ball $B_{n-1}$ in $\mathbb R^{n-1}.$ We calculate
$$\tag 1\frac{\partial \Sigma(y)}{\partial y_k} = (e_k,-y_k(1-|y|^2)^{-1/2}), \,\, k = 1,\dots , n-1.$$
Here $e_k$ is the standard basis vector in $\mathbb R^{n-1}.$
We want to think of the $n\times n$ matrix with $(1)$ giving the first $n-1$ column vectors, and $n(y),$ the unit normal vector, as the last column. But $n(y)$ is not $y,$ it is $\Sigma (y).$ ($\Sigma (y)$ is on the sphere; $y$ lives in $B_{n-1}.$) The determinant of this matrix is
$$\tag 2\frac{1}{(1-|y|^2)^{1/2}}.$$
Thus $d\sigma(\Sigma (y)) = dy/(1-|y|^2)^{1/2}.$ This is true in all dimensions, nicely enough.
To verify $(2),$ let's take advantage of the fact that $\Sigma$ is a graph. To review: Suppose $U\subset \mathbb R^{n-1}$ is open and $f:U\to \mathbb R$ is smooth. Then the graph $\{(y,f(y)):y\in U\}$ has surface area measure given by
$$d \sigma(y,f(y))= \left[(1+(D_1f(y))^2 + \cdots + (D_{n-1}f(y))^2\right]^{1/2}\,dy.$$
This follows from the general surface area formula you cited in your previous post. You should try to verify this. In the case of the upper hemisphere, $f(y) = (1-|y|^2)^{1/2}.$ A straightforward computation then gives $(2).$
Hopefully this will help with the trouble you are having with your specific integrals. Give it a try. I'll stop here for now.

Added later, another approach: Your integrals can be easily done if you accept one property of the surface measure $\sigma$ on the sphere: Rotation invariance. More precisely, if $T$ is an orthogal transformation of $\mathbb R^n$ and, say, $f$ is continuous on the unit sphere $S,$ then
$$\int_S f(y)\,d\sigma(y) = \int_S f(T(y))\,d\sigma (y).$$
So assume $i\in \{1,2\dots, n\},$ and $T$ is the orthogonal transformation that sends each standard basis vector $e_j$ to itself except for $j=i,$ where $T(e_i) = -e_i.$ Suppose $j\ne i,$ and $f(y) = y_iy_j.$ Then $f(T(y))=-f(y)$ for $y\in S.$ Thus
$$ \int_S y_iy_j\,d\sigma (y) = \int_S f(y)\,d\sigma(y) = \int_S f(T(y))\,d\sigma f(y) = -\int_S f(y)\,d\sigma(y).$$
Thus this integral is $0.$
For the other integral, suppose $i\ne j.$ Define $T$ to be the orthogonal transformation that switches $e_i$ with $e_j,$ and leaves the other basis vectors alone. Let $f(y)= y_i^2.$ Then
$$\int_S y_i^2\,d\sigma (y) = \int_S f(y)\,d\sigma (y) = \int_S f(T(y))\,d\sigma(y) =  \int_S y_j^2\,d\sigma (y).$$
From this we get, for fixed $i,$
$$\int_S y_i^2d\sigma (y) = \frac{1}{n}\int_S \left (\sum_{j=1}^{n}y_j^2\right )\,d\sigma (y) = \frac{1}{n}\int_S 1\,d\sigma (y) = \frac{\sigma (S)}{n}.$$
A: Solution Employing Symmetry
Let $S^{n-1}$ denote the unit hypersphere in $\mathbb{R}^n$ (with coordinate vector $\mathbf{y}=(y_1,y_2,\ldots,y_n)$) centered at the origin $\boldsymbol{0}_n$.  Write $\sigma_{n-1}$ for the hypersurface area measure of $S^{n-1}$.  Declare the northern hemisphere $S^{n-1}_+$ to be the one with $y_j>0$, and the southern hemisphere $S^{n-1}_-$ to be the one with $y_j<0$.  If $j\neq k$, then by symmetry, we have
$$\int_{S^{n-1}_-}\,y_jy_k\,\text{d}\sigma_{n-1}(\mathbf{y}) = -\int_{S^{n-1}_-}\,y_jy_k\,\text{d}\sigma_{n-1}(\mathbf{y})\,.$$
This shows that $$\int_{S^{n-1}}\,y_jy_k\,\text{d}\sigma_{n-1}(\mathbf{y})=\int_{S^{n-1}_+}\,y_jy_k\,\text{d}\sigma_{n-1}(\mathbf{y})+\int_{S^{n-1}_-}\,y_jy_k\,\text{d}\sigma_{n-1}(\mathbf{y})=0\,.$$
If $j=k$, then we note from symmetry that
$$\int_{S^{n-1}}\,y_1^2\,\text{d}\sigma_{n-1}(\mathbf{y})=\int_{S^{n-1}}\,y_2^2\,\text{d}\sigma_{n-1}(\mathbf{y})=\ldots=\int_{S^{n-1}}\,y_n^2\,\text{d}\sigma_{n-1}(\mathbf{y})\,.$$
Since $\sum\limits_{i=1}^n\,y_i^2=1$ on $S^{n-1}$,
$$\int_{S^{n-1}}\,y_j^2\,\text{d}\sigma_{n-1}(\mathbf{y})=\frac{1}{n}\,\int_{S^{n-1}}\,\sum_{i=1}^n\,y_i^2\,\text{d}\sigma_{n-1}(\mathbf{y})=\frac{1}{n}\,\int_{S^{n-1}}\,\text{d}\sigma_{n-1}(\mathbf{y})=\frac{1}{n}\,\Sigma_{n-1}\,,$$
where $\Sigma_{n-1}:=\displaystyle \int_{S^{n-1}}\,\text{d}\sigma_{n-1}(\mathbf{y})=\dfrac{2\pi^{\frac{n}{2}}}{\Gamma\left(\frac{n}{2}\right)}$ is the hypersurface area of $S^{n-1}$.  Here, $\Gamma$ is the usual gamma function.  That is, 
$$\frac{1}{n}\,\Sigma_{n-1}=\frac{\pi^{\frac{n}{2}}}{\Gamma\left(\frac{n}{2}+1\right)}\,,$$
which equals the volume of a unit $n$-dimensional hypersphere.

Solution Implementing the Divergence Theorem
Let $\lambda_n$ be the Lebesgue measure on $\mathbb{R}^n$ and $\mathbf{e}_1,\mathbf{e}_2,\ldots,\mathbf{e}_n$ the standard basis vectors of $\mathbb{R}^n$.  The normal vector to the hypersurface $S^{n-1}$ is given by $$\mathbf{n}=y_1\mathbf{e}_1+y_2\mathbf{e}_2+\ldots+y_n\mathbf{e}_n\,.$$  Write $B^n_r(\mathbf{x})$ for the open ball of radius $r>0$ centered at $\mathbf{x}\in\mathbb{R}^n$.  Note from the Divergence Theorem that $$\int_{\partial B^n_1(\boldsymbol{0}_n)}\,y_j\,y_k\,\text{d}\sigma_{n-1}(\mathbf{y})=\int_{\partial B^n_1(\boldsymbol{0}_n)}\,y_k\mathbf{e}_j\cdot\mathbf{n}\,\text{d}\sigma_{n-1}(\mathbf{y})=\int_{B^n_1(\boldsymbol{0}_n)}\,\big(\boldsymbol{\nabla}\cdot(y_k\mathbf{e}_j)\big)\,\text{d}\lambda_n(\mathbf{y})\,.$$   Obviously, $\boldsymbol{\nabla}\cdot(y_k\mathbf{e}_j)=0$ for $j\neq k$ and $\boldsymbol{\nabla}\cdot(y_j\mathbf{e}_j)=1$.
Consequently, the integral $\displaystyle\int_{S^{n-1}}\,y_jy_k\,\text{d}\sigma_{n-1}(\mathbf{y})$ equals $\lambda_n\big(B_1^n(\textbf{0}_n)\big)\,\delta_{j,k}$.  Here, $\delta$ is the Kronecker delta.  Since $\lambda_n\big(B^n_1(\textbf{0}_n)\big)=\dfrac{1}{n}\,\Sigma_{n-1}$, we get the same result as the first solution.
