Let me first describe where I start:
$$\iint_Sz^2\,dS$$
We want to compute the surface integral of the octant of a sphere $S$.
The radius = 1.
The sphere is centered at the origin.
$$S=x^2+y^2+z^2=1.$$
$$z=f(x,y)=\sqrt{1-x^2-y^2}$$ $R$ is the projection of $S$ on the $xy$-plane.
Now we compute the normalization factor used to project the integral on the $xy$-plane. First we compute the derivative $\frac{\partial f}{\partial x}$:
If
$$m=1-x^2-y^2$$
and
$$n=\sqrt{m}\,,$$
then the derivative of $n$ is
$$n'=\frac12m^{-\frac12},$$
and the derivative of $m$ (with regard to $x$) equal to:
$$m'=-2x.$$
Now we can compute $\frac{\partial f}{\partial x}$ using the chain rule:
$$\frac{\partial f}{\partial x}=n'\cdot m'=\frac12m^{-\frac12}\cdot m'=\frac12(1-x^2-y^2)^{-\frac12}\cdot-2x=-{\frac{x}{\sqrt{(1-x^2-y^2)}}}.$$
Since $\sqrt{(1-x^2-y^2)}=z$:
$$\frac{\partial f}{\partial x}=-{\frac xz}$$
And in the same manner (using the derivative of $m$ with regard to $y$) we can calculate $\frac{\partial f}{\partial y}$:
$$\frac{\partial f}{\partial y}=-{\frac yz}$$
Now we can calculate the normalization factor:
$$\sqrt{1+\left(\frac{\partial f}{\partial x}\right)^2+\left(\frac{\partial f}{\partial y}\right)^2}=\sqrt{1+\frac{x^2}{z^2}+\frac{y^2}{z^2}}=\frac1z\sqrt{x^2+y^2+z^2}.$$
Since $x^2+y^2+z^2=1$:
$$\sqrt{1+\left(\frac{\partial f}{\partial x}\right)^2+\left(\frac{\partial f}{\partial y}\right)^2}=\frac1z$$
Now we can calculate the projection $R$ of the sphere $S$ on the $xy$-plane: $$\iint_Sz^2\,dS=\iint_R z^2\frac1z \,dx\,dy=\iint_R z\;dx\,dy.$$ Substituting for $z$ finishes the conversion of the surface integral (remember $z=f(x,y)=\sqrt{1-x^2-y^2}$): $$\iint_Sz^2dS=\iint_R \sqrt{1-x^2-y^2}\;dx\,dy.$$
This is where my problem starts:
The book I'm reading says if we convert this to polar coordinates, the integration should be trivial.
So we convert to polar coordinates:
$$x=r\cos\theta$$
$$y=r\sin\theta$$
$$z=f(x,y)=f(r\cos\theta,r\sin\theta)$$
Calculate the Jacobian determinant:
$$\frac{\partial(x,y)}{\partial(r,\theta)}=\begin{vmatrix}\cos\theta & -r\sin\theta \\ \sin\theta & r\cos\theta\\ \end{vmatrix}=r\cos^2\theta+r\sin^2\theta=r\,(\cos^2\theta + \sin^2\theta)=r$$
Substitute:
$$\iint_Sz^2dS=\iint_R \sqrt{1-x^2-y^2}\;dx\,dy=\iint_T\sqrt{1-r^2\cos^2\theta-r^2\sin^2\theta}\cdot r \cdot dr\,d\theta$$
$$\iint_T\sqrt{1-r^2(\cos^2\theta+\sin^2\theta)}\cdot r \cdot dr\,d\theta = \iint_T\sqrt{1-r^2}\cdot r \cdot dr\,d\theta=\iint_T\sqrt{r^2-r^4}\cdot dr\,d\theta$$
Since the radius $r=1$ it is easy to see that:
$$\int^1_0\sqrt{r^2-r^4}\cdot dr=\int^1_0(r^2-r^4)^{\frac12}\cdot dr = \int^1_0 \frac{(r^2-r^4)^{\frac32}}{\frac32}=0$$
According to the book the result of the calculation of the surface of the sphere in the first octant should be $\pi/6$.
That won't happen if $\int^1_0\sqrt{r^2-r^4}\cdot dr=0$.
The domain of $\theta$ is:
$$0\le\theta\le\frac12\pi$$
So where am I going wrong?
It seems to me that you are using the power rule
$$\int (r^2 - r^4)^{1/2} dr= \frac{1}{3/2} (r^2 - r^4)^{3/2} +C,$$
but this is incorrect. Indeed if you differentiate the right hand side with respect to $r$, you do not get $(r^2 - r^4)^{1/2}$, but instead $(r^2 - r^4)^{1/2} (2r - 4r^3)$.
As noted in the comment, you can use substitution $u = 1-r^2$ to get
$$\int_0^1 \sqrt{1-r^2} rdr= -\frac 12 \int_1^0 \sqrt u du.$$
(we have $\int_1^0$ instead of $\int_0^1$ on the right hand side since you have to plug in the value accordingly: for $r=0$, you have $u= 1-r^2 = 1$, and for $r=1$ you get $u = 1-r^2 = 0$).
So you have
$$\int_0^1 \sqrt{1-r^2} dr = -\frac 12 \int_1^0 \sqrt u du = -\frac 13 u^{3/2}\bigg|_1^0 = \frac 13 $$
and thus $\int z^2 dS = \frac \pi 6$.