I have an isotropic 2D Harmonic Oscillator in cartesian coordinates \begin{equation} H = \frac{p_x^2}{2m} + \frac{p_y^2}{2m} + \frac{1}{2} m\omega^2 (x^2 + y^2) \end{equation} In terms of the usual creation and annihilation operators for the $x$ and $y$ modes, $n_x$ and $n_y$ this can be written as \begin{equation} H = \hbar\omega(n_x^\dagger n_x + n_y^\dagger n_y + 1) \end{equation} Now we can 'construct', three operators that commute with the Hamiltonian, apart from the (trivial) number operators $N_x$ and $N_y$: \begin{align} V_x &= a_x^\dagger a_y + a_y^\dagger a_x\\ V_y &= i(a_y^\dagger a_x - a_x^\dagger a_y)\\ V_z &= a_x^\dagger a_x - a_y^\dagger a_y \end{align} Now these operators can been proved to satisfy the commutation relations of angular momentum. In fact,
\begin{equation} V_i = a^\dagger \sigma_i a \end{equation} Where $a = \begin{pmatrix} a_x & a_y \end{pmatrix}$ and $\sigma_i$ are the Pauli matrices. This is really surprising to me as I don't see why and how this must be true. Can someone shed some light on this? I think for $2s+1$ dimensions, we can 'construct' such operators using the matrix representations of spin $s$. However, I haven't seen a proof of this.