Note: My approach lacks rigour and I'm rusty. I'm sure there are lots of holes, overly-complicated parts, and perhaps even just outright mistakes. However, I think the gist of the argument motivates why two operators which commute leads one to consider solutions of the form you mention.
Paraphrased: Why does does the fact that $\hat{A}$ and $\hat{B}$ commute in the Hamiltonian $\hat{H}=\hat{A} + \hat{B}$ mean that we can write the solution as factors of the solutions of $\hat{A}$ and $\hat{B}$?
Essentially, the fact that the Hamiltonian's parts, $\hat{A}$ and $\hat{B}$, commute means we can separate the general solution (of the form $\exp(-i\hat{H})$ into two exponential operators.
My approach lacks rigour, but I think it shows the idea:
For the equation
$\newcommand{\Ket}[1]{\left|#1\right>}$
$\newcommand{\Bra}[1]{\left<#1\right|}$
$\newcommand{\Braket}[2]{\left<{#1}|{#2}\right>}$
\begin{equation}
\hat{H}\Ket{\psi} = i\frac{\partial}{\partial t}\Ket{\psi} \tag{1}\label{SE}
\end{equation}
the general solution is known to be
\begin{align}
\Ket{\psi(t)} &= \exp\left(-it\hat{H} \right) \Ket{\psi(0)} \\
&= \exp\left(-it(\hat{A}+\hat{B}) \right) \Ket{\psi(0)}
\tag{2}
\end{align}
You can see this through substitution or look it up in any QM textbook or course notes.
The fact that $[A,B]=0$ means we can split the exponential operator, giving
\begin{align}
\Ket{\psi(t)} &= \exp\left(-it\hat{H} \right) \Ket{\psi(0)} \\
&= \exp\left(-it\hat{A} \right)\exp\left(-it\hat{B} \right) \Ket{\psi(0)}
\tag{3}
\end{align}
Let's try to write this in position representation.
First, let's recognise that we're dealing with more than one spatial coordinate and each operator is acting in mutually exclusive vector spaces, so that $\mathcal{H} = \mathcal{H}_A \otimes \mathcal{H}_B$. Alternatively, in position representation, we can see that $\hat{A}$ has an $x$-dependence and $\hat{B}$ has a $y$-dependence (and can be seen to commute by observation).
Also, the time-independent solutions for $\hat{A}$ and $\hat{B}$ separately in position representation are
\begin{align}
f(x) &= \Braket{x}{f} = \sum_i\lambda_i\Braket{x}{a_i} \\
g(y) &= \Braket{y}{g} = \sum_i\mu_i\Braket{x}{b_i}
\tag{4}
\end{align}
where the eigenfunctions are
\begin{align}
\hat{A}\Ket{a_i} &= a_i \Ket{a_i} \\
\hat{B}\Ket{b_i} &= b_i \Ket{b_i}
\tag{5}
\end{align}
Then using $\Braket{x,y}{\psi} = \psi(x,y)$ with $\Bra{x,y} = \Bra{x}\otimes\Bra{y}$, we get
\begin{align}
\Braket{x,y}{\psi} =& \left[\Bra{x}\exp\left(-it\hat{A} \right) \otimes \Bra{y}\exp\left(-it\hat{B} \right)\right] \Ket{\psi(0)}
\tag{6}
\end{align}
which we can write using the expansion of the identity as
\begin{align}
\Braket{x,y}{\psi} =& \left[\sum_n\Bra{x}\exp\left(-it\hat{A} \right) \Ket{a_n}\Bra{a_n} \right. \\
&\otimes \left. \sum_m\Bra{y}\exp\left(-it\hat{B} \right)\Ket{b_m}\Bra{b_m}\right]\psi(0)
\tag{7}
\end{align}
where the operators can now act on their respective eigenstates, giving
\begin{align}
\Braket{x,y}{\psi(t)} =& \sum_n\Braket{x}{a_n}\exp\left(-ita_n \right)\sum_m\Braket{y}{b_m}\exp\left(-itb_m \right) \\ &\int\int dx^\prime dy^\prime\Braket{a_n}{x^\prime}\Braket{b_m}{y^\prime} \Braket{x^\prime,y^\prime}{\psi(0)}
\tag{8}\label{FINALEQ}
\end{align}
where I've simplified the expression using $$\left(\Bra{x^\prime}\otimes\Bra{y^\prime}\right)\Ket{\psi(0)} = \Braket{x^\prime,y^\prime}{\psi(0)} \tag{9}$$
Note that $\Braket{x}{a_n}$ are the eigenfunctions of $\hat{A}$ and $\Braket{y}{a_m}$ are the eigenfunctions of $\hat{A}$ in position representation.
So, we can see that $\ref{FINALEQ}$ is simply a linear combination of factors of eigenstates of the individual operators $\hat{A}$ and $\hat{B}$:
\begin{align}
\Braket{x,y}{\psi(t)} =& \sum_{m,n}c_{m,n}f_n(x,t)g_m(y,t) \tag{10}
\end{align}
This is a general solution.
Ignoring the time-dependence, we can see that just one of those solutions (not constrained by any initial conditions) is $f_n(x,t)g_m(y,t) = \Braket{x}{a}\Braket{y}{b}$ i.e. the multiplication of the eigenstates of the operators separately. Generally, we can consider a linear combination of such and it's still a solution, and in this case, we've worked backwards from that.
I think this is also shows (maybe a bit strong given the lack of rigour... "suggests"?) that the independence of the variables leads us to consider states in different subspaces of $\mathcal{H}$, and hence operators in those different spaces commute, and therefore the general solution in position representation is formed of products of the separate solutions to those subspaces.
It doesn't matter that we're talking about $r$, $\phi$, and $\theta$, or something simpler like $x$ and $y$ - they're just specific representations of the fundamental states. Therefore, we can propose solutions of the form $\psi(x,y) = f(x)g(y)$.
Now, knowing this, whenever we see a Hamiltonian where $H=H_1 + H_2$, and the parts commute, we can jump straight to the conclusion and state $\psi(x,y) = f(x)g(y)$, where $f(x)$ and $g(y)$ are the solutions, and then immediately go about calculating $f(x)$ and $g(y)$.
Perhaps I'm misunderstanding but it seems the question is simply 'what are the eigenfunctions of a quantum particle confined to the (lateral) surface of a cylinder?'. There is no need to add an infinite potential to the Hamiltonian in order to confine the particle to the surface. For instance, when you solve the problem of a free particle in a 'box' in 2D, one does not add an infinite potential to confine the particle along the 3rd dimension. For simplicity, I will take the radius of the cylinder to be $R=1$. I'll also be assuming that question is about the scenario in which the cylinder is infinitely long.
There is no motion along the radial direction $r$, but there is motion along $z$ and $\phi$. The Hamiltonian given by Salmone is almost correct. It is ($\hbar=1$)
\begin{equation}
H = -\frac{1}{2m} \frac{\partial^2}{\partial z^2} -\frac{1}{2m} \frac{\partial^2}{\partial \phi^2},
\end{equation}
which is clearly separable. The eigenfunctions are then given by $\psi(z,\phi) = Z (z) \Phi(\phi)$, where $Z$ and $\Phi$ are the states of a free particle, i.e.,
\begin{align}
Z(z) &= e^{ikz}\\
\Phi(\phi) &= e^{ik\phi},
\end{align}
where I'll let you worry about the correct normalisation of these functions.
For $Z$, the values of $k$ are any real number, while for $\Phi$ we have $\Phi(\phi+2 \pi) = \Phi (\phi)$, which implies that $k = 0, \pm1,\pm2, \dots$.
Then the corresponding eigenvalues for either $Z$ and $\Phi$ are simply $E = k^2/2m$ (where $Z$ and $\Phi$ have different constraints on the possible values of $k$)
Hope this helps.
Best Answer
I think the other answer is (at least partially) mistaken. Let me first answer your question, and then explain why I have issues with the other answer.
Suppose you have a Hamiltonian of the form $H = H_1 + H_2$, where $[H_1, H_2] = 0$. Then, since $H_1$ and $H_2$ commute, they can be simultaneously diagonalized. That is, there exists an eigenbasis of the form $| \varepsilon_1, \varepsilon_2 \rangle$ where $H_1 | \varepsilon_1, \varepsilon_2 \rangle = \varepsilon_1 | \varepsilon_1, \varepsilon_2 \rangle$ and $H_2 | \varepsilon_1, \varepsilon_2 \rangle = \varepsilon_2 |\varepsilon_1, \varepsilon_2 \rangle$. An arbitrary state in your Hilbert space can be written in the form $$ | \psi \rangle = \sum_{\varepsilon_1, \varepsilon_2} c_{\varepsilon_1,\varepsilon_2} |\varepsilon_1, \varepsilon_2 \rangle , $$ but if your goal is only to find energy eigenstates, you are allowed to assume that you are working in a common eigenstate of $H_1$ and $H_2$ separately. This is the complete answer to your question. Note that if $H_1$ and $H_2$ didn't commute, then there would not exist a common eigenbasis, and it would not be possible to factorize your wavefunction.
Now, here is the problem with the other answer: there is no reason to assume that $H_1$ and $H_2$ have tensor product structure. Separation of variables works even when the Hamiltonian is NOT of the form $H = H_1 \otimes 1 + 1 \otimes H_2$. The difference is that when $H$ has this tensor product structure, the eigenvalues $\varepsilon_1$ and $\varepsilon_2$ can be chosen independently, but this doesn't have to be the case; indeed, this is not the case for the hydrogen atom!! In the hydrogen atom, the allowed quantum numbers $(n,\ell,m)$ cannot all be chosen independently, for example $\ell$ is not allowed to exceed $n$. This is indicative of the fact that, while the Hamiltonian splits into two commuting pieces, it does not have the tensor product structure indicated by the other answer.