Some details about the Box-Muller transform method I am confused about why $Z_0$ and $Z_1$ are independent. It seems like they both rely on $U_1$ and $U_2$. Could someone prove the statement?

 A: Let's not make any early conclusions about dependence or independence but rather try to obtain the transformed densities using basic principles. First, determine the range of the transformations. 
It is easy to see that the transformations are one-to-one and the set $$\left\{\left(u_1,u_2 \right): 0<u_1<1, 0<u_2<1 \right\}$$ is mapped onto 
$$\left\{\left(z_o, z_1 \right): -\infty< z_0 < \infty , -\infty < z_1 < \infty \right\}$$
Next we obtain the inverse transformations, that is solve for our original variables: $U_1$ and $U_2$. Divide $Z_1$ by $Z_0$ to obtain
$$\frac{Z_1}{Z_0} = \tan \left( 2\pi U_2 \right) \implies U_2 = \frac{1}{2\pi} \tan^{-1} \left( \frac{Z_1}{Z_0} \right) $$ 
Plugging that back in the first equation and using the fact that $\cos \left( \tan^{-1} (x) \right) = \frac{1}{\sqrt{1+x^2}}$ we can solve for $U_1$ to obtain
$$U_1 = \exp \left\{ -\frac{1}{2} \left( z_0^2 + z_1^2 \right) \right\}$$
It's of course good if you remember the properties of inverse trigonometric functions but I just saw it in wikipedia.
These transformations have the Jacobian
$$ J= \begin{bmatrix} - z_0 \exp \left\{ -\frac{1}{2} \left( z_0^2 + z_1^2 \right) \right\} & - z_1 \exp \left\{ -\frac{1}{2} \left( z_0^2 + z_1^2 \right) \right\} \\ - \frac{1}{2 \pi} \frac{z_1}{z_0^2 + z_1^2} & \frac{1}{2 \pi} \frac{z_0}{z_0^2 + z_1^2} \end{bmatrix}$$
the absolute determinant of which is simply 
$$|J| = \frac{1}{2 \pi} \exp \left\{ -\frac{1}{2} \left( z_0^2 + z_1^2 \right) \right\}$$
Recalling now that the joint pdf of two independent uniforms is the unit square, we have that the joint distribution of the transformatons is
$$f \left(z_0, z_1 \right) = \frac{1}{2 \pi } \exp \left\{ -\frac{1}{2} \left( z_0^2 + z_1^2 \right) \right\}, \ -\infty < z_0, \ z_1 <\infty$$
which we recognize as the the joint density of two independent standard normal variables. Nice, right?

It seems to me then that you forgot the Jacobian of the transformation, which makes the resulting pdf very non-constant. As a general remark, functional dependence does not necessarily imply statistical dependence. This is a recurring theme in statistics so be careful!
