Particle in a cylindrical well: role of $\varphi$? Look at a particle in a cylindrical well of radius $R$ and length $L$ with zero potential $V=0$ inside and $V=+\infty$ on the boundaries and outside the well. The time-independent Schrödinger Equation can be written as:
$$\nabla^{2} \psi=-k^{2}\psi\tag{1}$$
where $k^{2}=\frac{2mE}{\hbar^2}$
The geometry demands cylindrical coordinates, which yields from (1):
$$\frac{\psi_{\rho}}{\rho}+\psi_{\rho\rho}+\frac{\psi_{\varphi \varphi}}{\rho^{2}}+\psi_{zz}=-k^{2}\psi$$
for a wave function $\psi(\rho,\varphi,z)$, with boundary conditions:
$$\psi(\rho,\varphi,0)=0\\
\psi(\rho,\varphi,L)=0\\
\psi(R,\varphi,z)=0$$
Using the separation of variables Ansatz:
$$\psi(\rho,\varphi,z)=R(\rho)P(\varphi)Z(z),$$
this solves quite easily. For $P(\varphi)$ and noting that $\varphi$ is an angle, I get:
$$P(\varphi)=c_{1}\cos m\varphi$$
with $m=0,1,2,3...$.
$c_{1}$ is a normalisation constant that works out as:
$$c_1=\frac{1}{\sqrt{\pi}},$$
except for $m=0$, then:
$$c_1=\frac{1}{\sqrt{2\pi}}$$
For $Z$ and $R$ I get, respectively, a sine function and a Bessel of the first kind $J$.
So this is quite uneventful, really.

Then I tried a little 'check my work' and found this page on wave functions in cylindrical wells:
https://arxiv.org/abs/1205.3444
Slightly surprising, this page states the TISE as:
$$\frac{\psi_{\rho}}{\rho}+\psi_{\rho\rho}+\psi_{zz}=-k^{2}\psi,$$
with no $\varphi$ in sight. My first thought was, "Ah, symmetry!" but the angle $\theta$ does play in the wave functions of hydrogenic atoms and their potentials are symmetrical too.
Then I thought the presence of a radial potential may necessitate the angle, but the page I linked to also treats the case for a potential $V=\frac{\alpha}{r}$ (scroll down a little) and still no $\varphi$ mentioned.
I also solved the case of a circular potential well a while back and had to include the angle.
So my question is: does in the cylindrical case one have to include $P(\varphi)$ in the wave functions? And if not, why not?
 A: That paper mentioned in the question (arXiv:1205.3444) is, frankly speaking, very low quality.  It appears to be a conference contribution from a student who was working on programming Maple to solve the Schrödinger equation analytically.  However, there are two issues with the work.  The first, which you have identified, is that the paper misses many solutions.  The second, of less importance but still worth mentioning, is that the problem of solving the kind of boundary value problems analytically through a language like Maple was already solved at least a decade before this paper was posted on the arXiv.
So you are correct that the paper misses the dependence of the wave function on the azimuthal angle $\phi$.  There should indeed be an angular function, although your analysis is not quite complete in this respect either.  Using separation of variables, you should actually find for each integer $m\geq 0$, a general angular solution
$$P(\phi)=c_{1}\cos m\phi+c_{2}\sin m\phi.$$
Actually, it is usually more convenient to consider all integers $m$, positive or negative, with
$$P(\phi)=\frac{1}{\sqrt{2\pi}}e^{im\phi}.$$
This form is more convenient, because no special care needs to be taken for $m=0$; and moreover, in this formulation $m\hbar$ is the eigenvalue of the angular momentum component $L_{z}$.
Without the end cap boundary conditions at $z=0$ and $z=L$, the $z$-dependence of the wave function could be expressed either with trigonometric functions
$$Z(z)=d_{1}\cos(k_{3}z)+d_{2}\sin(k_{3}z)$$
(for $k_{3}\geq 0$), or complex exponentials
$$Z(z)=e^{ik_{3}z}$$
where any real $k_{3}$ is allowed and $\hbar k_{3}$ is the eigenvalue of the momentum component $p_{z}$. [In this case, the normalization of $Z(z)$ depends on what convention is being used to normalize continuum states, so I have not attempted to normalize the state here.]
With the end caps, only the
$$\sqrt{\frac{2}{L}}\sin(k_{3}z)$$
wave function is acceptable, and this makes $k$ subject to a quantization condition (just like $m$). The solution must have $k_{3}=\frac{n\pi}{L}$, and the wave functions are no longer eigenstates of $p_{z}$.
The separated equations for $R(\rho)$, $P(\phi)$ and $Z(z)$ are all second order; this is why there are two solutions for $P$ and $Z$ for each $m^{2}$ or $k_{3}^{2}$, respectively.  There should also be two solutions to the radial equation you found for $R(\rho)$; the two solutions are proportional to $J_{m}(k_{1}\rho)$ and $N_{m}(k_{1}\rho)$—where $k_{1}^{2}+k_{3}^{2}=k^{2}$—in terms of Bessel Functions of the first and second kind.  However, the Bessel Function of the second kind (or Neumann function) is not normalizable in the vicinity of the origin, so it does not get an acceptable solution. This leaves you with just the $J_{m}$ solution that you found.
Since $\rho=R$ must be a root of the Bessel Function, $k_{1}$ is quantized as well. $k_{1}R$ has to be one of the $z_{m,j}$, where $z_{m,j}$ is the $j$th root $J_{m}(z_{m,j})=0$ of the Bessel Function $J_{m}$. This makes the quantized energy level of the system
$$E_{j,m,n}=\frac{\hbar^{2}}{2m}\left(\frac{z_{m,j}^{2}}{R^{2}}+\frac{n^{2}\pi^{2}}{L^{2}}\right),$$
which depends on all three quantum numbers $j$, $m$, and $n$.
