Deriving expression for energy of a damped harmonic (1D) oscillator The equation
$$\ddot{x} + \gamma \dot{x} + \omega_0^2x=0 , \gamma \ll \omega_0$$
has the solution
$$e^{-(\gamma t)/2}(\alpha sin\ \omega_0t + \beta cos\omega_0t) = Ae^{-(\gamma t)/2} sin(\omega_0t + \phi)$$
(approximating $\frac{\sqrt{\gamma^2-4\omega_0^2}}{2}$ as $\omega_0i$ where $i$ is the imaginary unit)
The energy of this system is the sum of its kinetic and potential energy:
$$E(t)_ = \frac{1}{2}m\dot{x}^2 + \frac{1}{2}m\omega_0^2x^2$$
(using the product rule for $\dot{x}$ and then squaring that)
$$ = \frac{1}{2}mA^2e^{-\gamma t}\left[ \left( \omega_0cos(\lambda) -\frac{\gamma}{2}sin(\lambda)\right)^2+  \omega_0^2sin^2(\lambda)\right]$$
$$=\frac{1}{2}mA^2e^{-\gamma t}\left[\omega_0^2 +\left(\frac{\gamma^2}{4}sin^2(\lambda)-\omega_0 \gamma cos(\lambda) sin(\lambda)\right)\right]$$
where $\lambda$ = $\omega_0t + \phi$
I'm looking to end up with
$$\frac{1}{2}mA^2e^{-\gamma t}\omega_0^2$$
which is what is quoted in my lectures notes and from tutors, in order to attain later results about Q-factors equalling $\frac{\omega_0}{\gamma}$ et cetera.
But I can't see how I would get rid of the terms with $\gamma$ in the coefficient, unless I approximate both $\frac{\gamma^2}{4}$ and $\omega_0\gamma$ as $0$, which seems unrigorous. Any pointers would be much appreciated.
 A: Your answer is fine, even though I would not advise making such "unrigorous" approximations until you are a little more experienced, frankly. It's a little tedious (but not hard) to actually work out this problem completely and only make the approximation at the end. That gives the right answer too, and is a little easier to justify.
At any rate, this method did work this time. What you need to keep in mind is that you aren't just allowed to "throw away small numbers". The numbers that you throw away should be dimensionless for it to mean anything. ($c$ is a tiny number, if measured in parsecs per second.) The approximation tells you that $$\frac{\gamma}{\omega_0}\ll1,$$ so it's ok to ignore any number that is equal to or smaller than $\gamma/\omega_0$. Rearranging your own equation:
$$E = \frac{1}{2} m A^2 e^{-\gamma t} \omega_0^2 \left( 1 + \frac{\gamma^2}{4 \omega_0^2} \sin^2\lambda - \frac{\gamma}{2 \omega_0} \sin \left(2\lambda\right)  \right).$$
You can now ignore both the terms in the parentheses that contain $\gamma$, as your approximation says that they are very small. The result turns out to be what you'd expect:
$$E = \frac{1}{2} m A^2 e^{-\gamma t}\omega_0^2.$$
(To address the last line of your question, as well as your comment on the original post, it's not that you're throwing away $\omega_0 \gamma$, but rather that -- in the units chosen -- $\omega_0^2$ is so much larger than $\gamma^2$ and $\omega_0 \gamma$ that you can effectively ignore both these terms in comparison to $\omega_0^2$.)
