I have the Hamiltonian of an harmonic oscillator (with $\hbar=1$) $$ H = \omega \left(a^\dagger a + \dfrac{1}{2} \right) \;, $$ and the associated (canonical) partition function $$ Z = \text{Tr}\left[e^{-\beta H} \right] = \dfrac{1}{2} \text{csch}\left(\dfrac{\beta \omega}{2}\right) \;, $$ where $\beta = 1/T$ is the inverse temperature (with $k_B=1$).
For the operator $$ O = A(a+a^\dagger)+B(a+a^\dagger)^2 \;, $$ where $A$ and $B$ are real constants, I want to know its expectation value when evaluated on the harmonic oscillator at finite $\beta$.
Since this expectation value is given by $$ \langle O \rangle = \dfrac{\text{Tr}\left[O e^{-\beta H} \right]}{Z} \;, $$ the sum at the numerator can be easily computed if $O$ and $H$ can be diagonalised in the same basis. However, what should I do when this is not the case?