Solution (source code)

= Solution

Applying the <angular momentum operator> directly gives
$$
\hat L_z\psi=\hbar xz+i\hbar y f'(x)z.
$$
For an <eigenvalue> $\lambda$, comparison with $\lambda(f(x)+iy)z$ gives $\hbar x=\lambda f(x)$ and $\hbar f'(x)=\lambda$. The first identity excludes $\lambda=0$; differentiating it and using the second gives $\lambda^2=\hbar^2$. No additive constant in $f$ is permitted. \b[The possibilities are $f(x)=x$, $\lambda=\hbar$, and $f(x)=-x$, $\lambda=-\hbar$.]

For $f=x$, the polynomial $\psi=(x+iy)z$ is homogeneous of degree two and satisfies $\Delta\psi=0$. Expanding the three squared <orbital angular momentum> operators gives, with $D=x\partial_x+y\partial_y+z\partial_z$,
$$
\hat L^2=-\hbar^2\left(r^2\Delta-D^2-D\right).
$$
Here $D\psi=2\psi$ and $D^2\psi=4\psi$. Therefore \b[$\hat L^2\psi=6\hbar^2\psi$]. This is the degree-two case of <angular momentum of a homogeneous harmonic polynomial>; its angular labels are $l=2$, $m=1$.