A little prerequisite knowledge first, the author defines a 'mollifying' function $\rho_k: \mathbb{R}^n \rightarrow \mathbb{R}$ to have several properties, but that aren't relevant to my question (from what I understand, it exists to smooth solutions to PDEs, and enables us to show that the solution itself is smooth by definition). In our case, $\rho_k(x) = \rho_k(|x|)$ is radial.
Anyway, the author defines a harmonic function $u:\Omega \subset \mathbb{R}^n \rightarrow \mathbb{R}$ and then takes the 'convolution' of $u$ and $\rho_k$ to be: $$u_k(x) := (u \star \rho_k)(x) = \int_{\mathbb{R}^n}u(y)\rho_k(x-y)dy$$ After justifying that $u_k(x)$ is $C^{\infty}$, he goes on to show that $u_k = u$, and thus to deduce that $u$ is also $C^{\infty}$ and this is where I get lost; first he writes: $$u_k(x) = \int_{\mathbb{R}^n}u(y)\rho_k(x-y)dy = \int_{\mathbb{R}^n}u(x-z)\rho_k(z)dz = \int_{\mathbb{R}^n}u(x-z)\rho_k(|z|)dz$$ Which I understand; however then, he states that we need to transform $z$ into polar coordinates, from which he subsequently writes: $$\begin{align}\int_{\mathbb{R}^n}u(x-z)\rho_k(|z|)dz = \int_{0}^{\infty}\int_{S^{n-1}}u(x + z(r,\theta))\rho_k(r)r^{n-1}d\theta dr &= \int_{S^{n-1}}\int_{0}^{\infty}u(x)\rho_k(r)r^{n-1}d\theta dr\end{align}$$ First, what substitution exactly is made here? And what is $S^{n-1}$? I presume it is the unit sphere in $\mathbb{R}^{n-1}$, but then why doesn't he just write the limits as $(0,2\pi)$? Also, why is an initial integral over $\mathbb{R}^n$ equivalent to taking this kind of product integral? Also I think he uses the mean value property to get from the second to third bit, but I'm not exactly sure how.
It's all confusing me a little and I just want a bit of clarification for, what the author would probably tell me are, 'obvious' points...
