Rationale for writing wave function as product of independent wave functions

It's because when $$V(x,y,z) = V_x(x) + V_y(y) + V_z(z),$$ (I guess that your extra identity $V(x,y,z)=V(z)$ is a mistake), we also have $$ H = H_x + H_y + H_z$$ because $H = (\vec p)^2 / 2m + V(x,y,z) $ and $(\vec p)^2 = p_x^2+p_y^2+p_z^2$ decomposes to three pieces as well.

One may also see that the terms such as $H_x\equiv p_x^2/2m+V_x(x)$ commute with each other, $$ [H_x,H_y]=0 $$ and similarly for the $xz$ and $yz$ pairs. That's because the commutators are only nonzero if we consider positions and momenta in the same direction ($x$, $y$, or $z$).

At the end, we want to look for the eigenstates of the Hamiltonian $$ H|\psi\rangle = E |\psi \rangle$$ and because we have $H = H_x+H_y+H_z$, a Hamiltonian composed of three commuting pieces, we may simultaneously diagonalize them i.e. look for the common eigenstates of $H_x,H_y,H_z$, and therefore also $H$. So given the separation condition for the potential, we may also assume $$ H_x |\psi\rangle = E_x |\psi\rangle $$ and similarly for the $y,z$ components. However, the equation above is just a 1-dimensional problem that implies that $|\psi\rangle$ must depend on $x$ as a one-dimensional quantum mechanical energy eigenstate wave function, $$ \psi(x) = C\cdot \psi_n(x) $$ which is an eigenstate of $H_x$. This has to hold but the normalization factor is undetermined. We usually say that it's a constant but this statement only means that it is independent of $x$. In reality, it may depend on all observables that are not $x$ such as $y,z$. So a more accurate implication of the $H_x$ eigenstate equation is $$ \psi(x,y,z) = C_x(y,z)\cdot \psi_{n_x}(x) .$$ In a similar way, we may show that $$ \psi(x,y,z) = C_y(x,z)\cdot \psi_{n_y}(y) $$ and $$ \psi(x,y,z) = C_z(x,y)\cdot \psi_{n_z}(z) $$ and by combining these three formulae, we see that the whole function must factorize to a product of functions of $x$ and $y$ and $z$ separately. If you need a rigorous proof of the last simple step, take e.g. the complex logarithms of the three forms for $\psi$ above and compare e.g. the first pair: $$\ln\psi = \ln C_x(y,z) +\ln\psi_{n_x}(x) = \ln C_y(x,z)+\ln \psi_{n_y}(y) $$ Take e.g. the partial derivative of the last equation with respect to $y$: $$ \frac{\partial \ln C_x(y,z)}{\partial y} = \frac{\partial \ln\psi_{n_y}(y) }{ \partial y }$$ The other two (1+1) terms are zero because they didn't depend on $y$. The right hand side above only depends on $y$, so the same must be true for the left hand side. I am going to make a simple conclusion but to make it really transparent, let's differentiate the latter equation over $z$, too. The $\psi_{n_y}$ term disappears as well so we have $$\frac{\partial^2 \ln C_x(y,z)}{\partial y\,\partial z} = 0$$ It means that $\ln C_x(y,z)$ must have the form $K_x(y)+L_x(z)$, and $e^{K_x(y)}e^{L_x(z)}$ must be the remaining factors in the wave function.

We say that the wave function in the product form is a "tensor product" of the three independent one-dimensional wave functions and more "operationally", as another user mentioned, the method described above is the method of "separation of variables".

This method is called "separation of variables", and it is one of several strategies for finding solutions to multi-dimensional field problems in physics.

It's big advantage is that solving N 1-dimensional differential equations is generally easier than solving 1 N-dimensional problems (N>1), but it is contingent on there being a "uniqueness theorem" for the category of problems that you are looking at. Happily many common field problems in physics have such a theorem.

To show that the condition on the potential is required for the Schrödinger equation simply make the separation, and expand. If the potential has the above form, you can write the LHS of the equation in three terms: $$ \left[\left( \frac{-\hbar^2}{2m} \frac{\partial^2}{\partial x} + V(x) \right)\psi(x)\right] \left[\left( \frac{-\hbar^2}{2m} \frac{\partial^2}{\partial y} + V(y) \right)\psi(y)\right] \left[\left( \frac{-\hbar^2}{2m} \frac{\partial^2}{\partial z} + V(z) \right)\psi(z)\right] $$ and you can clearly rewrite the whole into three separate conditions like $$ \left( \frac{-\hbar^2}{2m} \frac{\partial^2}{\partial x} + V(x) \right)\psi(x) = i\hbar{}\frac{\partial}{\partial t} \psi(x) . $$ On the other had, if the potential can not be written in this way, you can not get the LHS into the form you need and you can not proceed along these lines.