Let's quickly review the quantum harmonic oscillator. We have a single particle moving in one dimension, so the Hilbert space is $L^2(\mathbb{R})$: the set of square-integrable complex functions on $\mathbb{R}$. The harmonic oscillator Hamiltonian is given by
$$H= \frac{P^2}{2m} + \frac{m\omega^2}{2}X^2$$
where $X$ and $P$ are the usual position and momentum operators: acting on a wavefunction $\psi(x)$ they are $X \psi(x) = x\psi(x)$ and $P \psi(x) = -i\hbar\ \partial \psi / \partial x$. Of course, we can also think of them as acting on an abstract vector $|\psi\rangle$.
By letting $P \to -i\hbar\ \partial/\partial x$ we could solve the time independent Schrödinger equation $H \psi = E \psi$, but this is a bit of a drag. So instead we define operators $a$ and $a^\dagger$ as in your post. Notice that the definition of $a$ and $a^\dagger$ has nothing whatsoever to do with our Hamiltonian. It just so happen that these definitions are convenient because the Hamiltonian turns out to be $\hbar \omega (a^\dagger a + 1/2)$.
For convenience we define the number operator $N = a^\dagger a$; at this stage number is just a name with no physical interpretation. Using the commutation relation $[a,a^\dagger] = 1$ and some algebra we notice that $N$ has a nondegenerate spectrum given by the natural numbers. In other words, the eigenvalues of $N$ are $\{0,1,2,\dots\}$, and to each eigenvalue $n$ there corresponds a single state $|n\rangle$ with $N|n\rangle = n |n\rangle$. Notice that, again, $N$ is independent of our Hamiltonian. However, because the Hamiltonian turns out to be $\hbar \omega (N+1/2)$ we immediately know that the states $|n\rangle$ are its eigenvectors, with energies $\hbar \omega (n + 1/2)$.
Now you are given a different Hamiltonian. The Hilbert space is still exactly the same, and so are $a$, $a^\dagger$ and $N$, because their definition had nothing to do with the original Hamiltonian. You can still use their properties to find energies, eigenvectors, and so on. The states $|n\rangle$ are still the eigenstates of $N$, though a priori they might not be eigenstates of the new $H$ (exercise 31 asks you to prove that they in fact are eigenstates of the new $H$). The important point here is that operators are (usually) defined independently of the Hamiltonian. They characterize the physical system. After all, you know that there are operators $X$ and $P$, and you have no qualms about using them with different Hamiltonians. The Hamiltonian gives the energy and the time evolution, but the observables and related operators are independent of your choice of Hamiltonian.
About the physical interpretation... exercise 31 asks you to prove that $H=\hbar\omega_0 N + \hbar \omega_1 (N^2-N)$; notice that we have gotten rid of $\hbar\omega_0 /2$ since it is just a constant. I would usually expect $\omega_1$ to be smaller than $\omega_0$ so this is a small perturbation (for small $n$ at least), but we don't really care about that right now. You can see that $|n\rangle$ are still the eigenstates of the Hamiltonian; all we did is shift the energies by an amount $\hbar \omega_1 (N^2-N)$.
First, you need to understand $a_i$ is a bosonic operator satisfying the bosonic commutation relation $[a_i,a_j]=0$, meaning that $a_ia_j=a_ja_i$. Now we show that $\epsilon_{ijk}a_ja_k=0$. Because if you fix $i=1$, then
$$\epsilon_{1jk}a_ja_k=\epsilon_{123}a_2a_3+\epsilon_{132}a_3a_2=a_2a_3-a_3a_2=0.$$
For other choice of $i$, the proof is similar. Then $\epsilon_{ijk}a_ja_k=0$ implies $\epsilon_{ijk}a_j^\dagger a_k^\dagger=0$ by Hermitian conjugating both sides of the equation. Because both $\epsilon_{ijk}a_ja_k$ and $\epsilon_{ijk}a_j^\dagger a_k^\dagger$ are zero, so their difference is also zero.
Best Answer
The "cartesian" basis for the 2D Harmonic oscillator is given in terms of the $x$ and $y$ eigenstates of the corresponding number operators $\hat n_x$ and $\hat n_y$, written as $|n_x\rangle \otimes|n_y\rangle$. These number operators have to been understood as $\hat n_x \equiv \hat n_x \otimes 1_y$ and $\hat n_y \equiv 1_y\otimes \hat n_y $, with $1_x$ and $1_y$ the identity matrix in the $x$ and $y$ space.
One can also solve the same problem using another basis, spanned by the $g$ and $d$ operators $$\hat a_d=\frac{1}{\sqrt{2}}\left(\hat a_x \otimes 1_y - i1_x \otimes \hat a_y\right),\\ \hat a_g=\frac{1}{\sqrt{2}}\left(\hat a_x \otimes 1_y + i1_x \otimes \hat a_y\right),$$ which act on both $x$ and $y$ spaces at the same time. The effect of the $g$ operator is $$\hat a_g |n_x\rangle \otimes|n_y\rangle = \frac{1}{\sqrt{2}}\left(\sqrt{n_x}|n_x-1\rangle \otimes|n_y\rangle+\sqrt{n_y}|n_x\rangle \otimes|n_y-1\rangle\right).$$
One can then show that the Hamiltonian can be expressed in terms of the $g$ and $d$ number operators (in particular $\hat n_x\otimes 1_y+1_x \otimes\hat n_y= \hat n_d + \hat n_g$).
In principle, one can express that eigenstates of $\hat n_d$ and $\hat n_g$, $|n_g,n_d\rangle$ in the basis $|n_x\rangle \otimes|n_y\rangle$.