Solution (source code)

= Solution

Every <orbital angular momentum> operator annihilates a radial function: substituting $\partial_i f(r)=f'(r)x_i/r$ makes the two terms cancel. It also annihilates $r^2$. Thus the radial factor can be carried through the angular operators. Since $x_3^2$ is independent of $x_1,x_2$, $L_3\psi=0$ for every value of $a$.

Direct differentiation gives
$$
L_1^2x_3^2=2\hbar^2(x_3^2-x_2^2),\qquad
L_2^2x_3^2=2\hbar^2(x_3^2-x_1^2),\qquad L_3^2x_3^2=0.
$$
Consequently $L^2x_3^2=\hbar^2(6x_3^2-2r^2)$ and $L^2r^2=0$. For a nonzero admissible radial factor, comparison with $\lambda(x_3^2+ar^2)f(r)$ forces $\lambda=6\hbar^2$ and $6a=-2$. Therefore \b[the simultaneous eigenfunction occurs at]
$$
\boxed{a=-\frac13,\qquad L^2\psi=6\hbar^2\psi,\qquad L_3\psi=0.}
$$
It is the <axially symmetric quadratic orbital eigenfunction>, with quantum numbers $l=2$, $m=0$ and angular dependence proportional to the <spherical harmonic> $Y_{20}$. The subtraction of $r^2/3$ removes its rotationally invariant scalar component.