Radial quantum number for infinite circular well For completeness, I will sketch the solution of a particle in an infinite circular well first and then get to my question. I apologize in advance since the introduction is standard undergraduate quantum mechanics.

Let us consider the problem of a particle in an infinite circular well. This is described by the potential
$$V(r) = \begin{cases}0&r \leq R \\ \infty&r> R \end{cases}$$
The Schrödinger equation then splits into a radial and an angular part and we can write the eigen-function as
$$\psi(r,\theta) = u(r) e^{i l \theta}$$
where $l = 0,\pm 1,\pm2,..$ due to the single-valuedness of the wave-function. The radial part of the equation is
$$\left(\frac{\partial^2}{\partial r^2} + \frac{1}{r}\frac{\partial}{\partial r} - \frac{l^2}{r^2} + \frac{2mE}{\hbar^2}\right)u(r) = 0$$
with the boundary condition: $u(R) = 0$. The solutions to this are given by the regular Bessel functions (we discard $Y_l$'s since they blow up at r = 0)
$$u_{n,l}(r) = J_{l}\left(\alpha_{n,l}\frac{r}{R}\right)$$ 
where $\alpha_{n,l}$ is the $n^{th}$ zero of $J_{l}(r)$.
In this way, we can construct normalizable wavefunctions
$$\psi_{n,l}(r,\theta) = \mathcal{N}_{n,l}J_{l}\left(\alpha_{n,l}\frac{r}{R}\right) e^{il\theta}$$
where $\mathcal{N}_{n,l} = \frac{1}{\sqrt{\pi}R |J_{l+1}(\alpha_{n,l})|}$. These are simultaneous eigenfunctions of the Hamiltonian and of angular momentum, with energy $E_{n,l} = \frac{\hbar^2 \alpha_{n,l}^2}{2mR^2}$ and angular momentum = $l \hbar$.

Now, my question is the following: Suppose I want to evaluate 
$$(\partial_x \pm i \partial_y) \psi_{n,l} = e^{\pm i \theta}\left(\frac{\partial}{\partial r} \pm \frac{i}{r}\frac{\partial}{\partial \theta} \right)\psi_{n,l}(r,\theta)$$
Using the properties of Bessel functions, I find that this evaluates to
$$\mp\mathcal{N}_{n,l}\frac{\alpha_{n,l}}{R}e^{i(l\pm 1)\theta}J_{l\pm 1}\left(\alpha_{n,l}\frac{r}{R} \right)$$
This is no longer an eigenfunction since $J_{l\pm 1}\left(\alpha_{n,l}\frac{r}{R} \right) \neq 0$ when $r = R$. 
So, is there any way of constructing eigenfunctions with fixed angular momentum, $l$, such that when I act on them by the operator(s) $\partial_x \pm i \partial_y$, I get another eigenfunction? 
 A: Define $D=\{(x,y)\in\mathbb R^2:\|(x,y)\|\leq R\}$ as the disk of interest. There are two spaces of interest here: the space of square-integrable functions on $D$, $L_2(D)$, and the space of such functions with Dirichlet boundary conditions, $\mathcal H=\{\psi\in L_2(D):\psi(p)=0\:\forall p\in \partial D\}$.
You're interested in the hamiltonian $H=\frac{1}{2m}\hat{\mathbf p}^2$, with domain 
$$\mathscr D(H)=\{\psi\in\mathcal H:H\psi\in L_2(D)\},$$
and the angular momentum operator $L=xp_y-yp_x$ with domain 
$$\mathscr D(L)=\{\psi\in\mathcal H:L\psi\in L_2(D)\}.$$
As you've shown, the only simultaneous eigenfunctions of these two operators are given by
$$\psi_{n,l}(r,\theta) = \mathcal{N}_{n,l}J_{l}\left(\alpha_{n,l}\frac{r}{R}\right) e^{il\theta}.$$
As such, there is no way to construct simultaneous eigenfunctions of $H$ and $L$ with a different behaviour for the prospective ladder operator - the theory is completely constrained.
The relationship between the eigenfunctions of $H$ and the momentum operators is complex and full of subtle details. There's nothing particularly new here; you'll get much the same issues with one-dimensional momentum in a finite 1D well. For more details see e.g. this paper.

If I understand your question correctly, you're looking for eigenfunctions of $L$ which are not necessarily eigenfunctions of $H$ such that $\partial_x \pm i \partial_y$ will act as a ladder operator. This is an easier problem to pose. You're then looking for a basis set of $\phi_{nl}\in \mathcal H$ of the form $\phi_{nl}(r,\theta)=f_{nl}(r) e^{il\theta}$ for which $(\partial_x\pm i\partial_y) \phi_{nl}=A_{nl}^\pm \phi_{n,l\pm1}$, or in other words
$$
\partial_l^\pm f_{nl}(r)=\left(\frac{\partial}{\partial r} \mp \frac{l}{r} \right)f_{nl}(r) 
=A_{nl}^\pm f_{n,l\pm 1}(r) 
\qquad\text{under}\ \ f_{nl}(R)=0.
$$
For this to happen you want to impose a consistency requirement - going one step up and then back again should land you in the same spot as before, i.e.
\begin{align}
\partial_{l+1}^- \partial_l^+ f_{nl}(r) 
 = \partial_{l+1}^-A_{nl}^+f_{n,l+1}(r) 
 = A_{n,l+1}^-A_{nl}^+f_{nl}(r) 
 = -B_{nl}f_{nl}(r) 
\end{align}
which reduces to the form
$$
\left(r^2\frac{\partial^2}{\partial r^2}+r\frac{\partial}{\partial r} + (B_{nl}r^2-l^2)\right)f_{nl}=0.
$$
This is Bessel's equation, which you'd expect because nothing so far uses the finite boundary so it needs to accommodate the case with no such boundary, for which the Bessel solutions work.
This consistency requirement has two curious features. One is that it doesn't matter which way you go, and indeed a short calculation shows that
$$\partial_{l+1}^- \partial_l^+=\partial_{l-1}^+ \partial_l^-.$$
The second follows from this, and it's that the consistency requirement is all you really need, and if $f_{nl}$ is a solution for one $l$ then you can use it to build other states $\partial_{l}^\pm f_{nl}$ and these will still obey the consistency requirement.
The problem is, though, that this is not compatible with the finite-disk boundary conditions. Indeed, you want the $\phi_{nl}$ to remain in $L_2(D)$, and this rules out any component along the Neumann-Weber solution $Y_l$, so $f_{nl}\propto J_l$. Asking for $\phi_{nl}\in \mathcal H$ then forces you back to the eigenfunctions of $H$, and you know that you can only do this for one eigenfunction at a time - you can impose it for $f_{nl}$ but it will then necessarily break for $\partial_l^\pm f_{nl}$.
So I guess the short answer is simply "No".
A: It is straight forward to show that it is not possible to find eigenfunctions of $\hat L_z$ that are also eigenfunctions of $\partial_x \pm i\partial_y$ as these operators do not commute.
Indeed, using $\hat L_z = \hat x \hat p_y - \hat y \hat p_x $ and $\partial_x \pm i\partial_y=-i \hat p_x \mp \hat p_y$, as well as the standard canonical commutation relations, one finds
$$
[\hat L_z,\partial_x \pm i\partial_y]=\mp(\partial_x \pm i\partial_y).
$$
