You have an isotropic three dimensional quantum harmonic oscillator so the Hamiltonian is $$ H=\frac{p^2}{2}+\frac{r^2}2 $$ If you do the creation-annihilation operator-algebra trick and define creation operator $$ a^\dagger_j=\frac{x_j-ip_j}{\sqrt{2\hbar}} $$ so that the Hamiltonian is $$ H=\hbar \left (a_1^\dagger a_1+a_2^\dagger a_2+a_3^\dagger a_3+\frac 3 2 \right) $$ so that the Hamiltonian is diagonalized by states characterized by the excitation number in each direction.
BUT if we want to compute angular momentum, we could define $3$-rotating creation operators $$ a_\uparrow^\dagger=\frac{a^\dagger_1+ia^\dagger_2}{\sqrt 2} \ \ \ \ \ \ \ \ \ a_\downarrow^\dagger=\frac{a^\dagger_1-ia^\dagger_2}{\sqrt 2} $$ and, as $a_1^\dagger a_1+a_2^\dagger a_2=a_\uparrow^\dagger a_\uparrow+a_\downarrow^\dagger a_\downarrow$, the Hamiltonian is still nicely diagonized; moreover we also have the 3-angular momentum operator diagonal, $L_3=\hbar (a_\uparrow^\dagger a_\uparrow-a_\downarrow^\dagger a_\downarrow)$, so $L_3 \left | n_\uparrow n_\downarrow n_z \right \rangle=\hbar (n_\uparrow-n_\downarrow) \left | n_\uparrow n_\downarrow n_z \right \rangle$.
How might you similarly diagonalize $L^2$ with respect to creation-annihilation? The spherical harmonics are hellish, but Schwinger's creation operator approach is so nice and intuitive. How could we compute $L^2$ probabilities without ever integrating again?
Edit More precisely: Are there commuting operators $a,b,c$ such that $[a,a^\dagger]=\hbar$ (etc for $b$ and $c$) and the set of states defined by $$\left | n_a n_b n_c \right \rangle=(a^\dagger)^{n_a}(b^\dagger)^{n_b}(c^\dagger)^{n_c}\left | \psi_0 \right \rangle$$ form an orthogonal basis diagonalizing $H,L^2$,and $L_3$? Here $\left | \psi_0 \right \rangle$ is the ground state. How can you determine $a,b,c$ algebraically?