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)$.
Congratulations! You found out that the time dependence of the harmonic oscillator's eigenstates do not resemble the classical oscillator. If you want a non-zero expectation value you should prepare the system in a superposition of adjacent eigenstates, like
$$
|\psi\rangle = c_0 |0\rangle + c_1|1\rangle.
$$
That's a consequence of $x$ depending on $a + a^\dagger$.
Either way, if you want the state that truely resembles the classical oscillator you should look at the coherent states. There are many ways to define them, one example that makes clear their resemblance to the classical oscillator is to translate by a finite distance $d$ the ground state:
$$
|\psi\rangle = \exp \left (-\frac{i p d}{\hbar} \right )|0\rangle.
$$
Using the Heisenberg picture, where the time-dependent operator $x$ is
$$
x(t) = x(0) \cos \omega t+\frac{p(0)}{m\omega} \sin \omega t
$$
and $|\psi\rangle$ is fixed on time,
you can prove that the expectation value of $x(t)$ evolves just like a classical oscillator of amplitude $d$:
$$
\langle x(t)\rangle = \langle \psi |x(t)|\psi\rangle = d \cos \omega t.
$$
Best Answer
The trace can be calculated in any basis. Therefore, we will calculate it in the oscillator eigenbasis $\{|n\rangle\}$ as follows. First, $$\text{Tr}\left( e^{-H/ k_B T} \right) = \sum_{n=0}^{\infty}\langle n|e^{-H/ k_B T}|n\rangle =\sum_{n=0}^{\infty}\langle n|e^{-\hbar\omega(n+1/2)/ k_B T}|n\rangle =\sum_{n=0}^{\infty}e^{-\hbar\omega(n+1/2)/ k_B T}. $$ This last sum can be evaluated because it can be written as a geometric series: $$\text{Tr}\left( e^{-H/ k_B T} \right) =e^{-\hbar\omega/2k_B T}\sum_{n=0}^{\infty}\left(e^{-\hbar\omega / k_B T}\right)^n =e^{-\hbar\omega/2k_B T}\frac{1}{1-e^{-\hbar\omega / k_B T}} =\frac{e^{\hbar\omega / 2k_B T}}{e^{\hbar\omega / k_B T}-1} $$ Then, note that $$\text{Tr}(\rho H) = \text{Tr}\left(\frac{e^{-H / k_B T}}{\text{Tr}\left( e^{-H/ k_B T} \right) } H\right) =\frac{\text{Tr}\left(e^{-H / k_B T}H\right)}{\text{Tr}\left( e^{-H/ k_B T} \right) }. $$ So, we have to do the top: \begin{align*} \text{Tr}\left(e^{-H / k_B T}H\right) &=\sum_{n=0}^{\infty}\langle n|e^{-H / k_B T}H|n\rangle =\sum_{n=0}^{\infty}e^{-\hbar\omega(n+1/2)/ k_B T}\hbar\omega\left(n+\frac{1}{2}\right)\\ &=\frac{\hbar\omega}{2}\sum_{n=0}^{\infty}e^{-\hbar\omega(n+1/2)/ k_B T} + \hbar\omega\sum_{n=0}^{\infty}e^{-\hbar\omega(n+1/2)/ k_B T}n\\ &=\frac{\hbar\omega}{2}\frac{e^{\hbar\omega / 2k_B T}}{e^{\hbar\omega / k_B T}-1} + \hbar\omega \frac{e^{\hbar\omega / 2k_B T}}{e^{\hbar\omega / k_B T}-1} \frac{1}{e^{\hbar\omega / k_B T}-1} \end{align*} Putting these together yields $$\text{Tr} (\rho H) = \frac{\hbar \omega}{2} + \frac{\hbar \omega}{e^{\hbar \omega / k_B T} - 1} \, .$$