Calculating the Variance of the Thermal state

We have a Harmonic Oscillator in the Thermal state $tau(beta)$ which is defined

$$tau(beta) = frac{e^{-beta H}}{mathrm{Tr}(e^{-beta H})}$$

where $Z = mathrm{Tr}(e^{-beta H})$ is known as the partition function.

Now I was asked to calculate the average initial energy

$$E(tau(beta)) = mathrm{Tr}(tau H)$$

Now I did this (correctly) by doing:

$$E(tau(beta)) = (frac{1}{Z}) sum_n n cdot mathrm{Exp}[-n beta]$$

Now I was asked to calculate the corresponding variance.
I am supposed to obtain

$$V(tau(beta)) =(frac{1}{Z}) sum_n (n -varepsilon)^2 cdot mathrm{Exp}[-n beta] tag{1}$$

where I am not 100% sure what $varepsilon$ is. I used the formula

$$V(rho) = mathrm{Tr}(rho H^2) – (mathrm{Tr}(rho H))^2$$

but am not able to understand how I am supposed to get eq. (1) with this. Can anyone help?