Quantum Harmonic Oscillator and the Classical Limit We can solve for the stationary states of a quantum harmonic oscillator denoted by $|n\rangle$ with energy eigenvalues $(n+\frac{1}2)\hbar\omega$. However if our system is in a stationary state, the time evolution is given by $|\psi,t\rangle=e^{\frac{-iE_nt}{\hbar}}|\psi,0\rangle$ which translates to $|\psi,t\rangle=e^{-i(n+\frac{1}2)\omega t}|n\rangle$ here. Then $\langle x \rangle=\langle\psi|x|\psi\rangle=0$ (because the time dependence drops out from the conjugation and we can show that $\langle x\rangle=0$ for all stationary states), and consequently we do not retrieve the classical sinusoidal motion.
However it can be shown that if we have $|\psi,t\rangle=\displaystyle\sum_{n}a_ne^{-(n+\frac{1}2)\omega t}|n\rangle$ such that more than one $a_n$ is non-negligible, we obtain $\langle x\rangle=\displaystyle\sum_{n}X_n\cos(\omega t+\phi_n)$ which is the classical sinusoidal motion we expect.
So I conclude that for classical motion we have to be in a superposition of stationary states. However, I have seen many arguments where the classical limit of the quantum harmonic oscillator is considered simply by looking at the form of the probability distribution for large $n$, for example here http://hyperphysics.phy-astr.gsu.edu/hbase/quantum/hosc6.html (check out the graphs towards the bottom). This seems to be contradicting my above analysis of the situation - if we have to be in a superposition of stationary states to retrieve classical motion, why bother considering the large $n$ limit of stationary states - surely this information is useless in the classical limit because the classical limit never involves simply one stationary state - if it did (as I worked out above) there would be no sinusoidal motion.
Thanks for any help :)
 A: First of all note that what you have found is perfectly consistent with classical mechanics (CM) because the expectation value of a classical oscillator is also $x=0$ because of the symmetry. Yes you have to have a (particular) superposition of energy eigenstates in order to get the classical motion. The Heisenberg representation of the QM is particularly useful for this task. For example the ground state translated by an amount $a$ behaves exactly like the classical oscillator. There is a nice applet by falstad that you can play and see the motion by yourself if you don't know anything about the coherent states or Heisenberg representation. If you know the Heisenberg representation then I guess you and I both agree that the following looks classical:
\begin{align}
\hat{x_H}(t) &= \hat{x} \cos(\omega t) + \frac{\hat{p}}{m \omega} \sin(\omega t)\\
\hat{p_H}(t) &= \hat{p} \cos(\omega t) - m\omega \hat{x} \sin(\omega t)
\end{align}
However you should also know that there is no easy transition from QM to CM. 
A: 
This seems to be contradicting my above analysis of the situation - if
  we have to be in a superposition of stationary states to retrieve
  classical motion, why bother considering the large n limit of
  stationary states

There is no contradiction. 
In the hyperphysics link you gave, they show that you need n>>1 to recover the classical probability. However that is not all you need. It is necessary, not sufficient.
Notice that it is not true that simply taking n>>1 gives you the classical probability. The n>>1 probability (as exemplified by $|\psi_{10}|^2(x)$ in the link you gave) is not actually the same (nor is it really even similar) to the classical probability. Rather, $|\psi_{10}|^2(x)$ oscillates up and down a lot, but has an envelope that resembles the classical probability... 
The authors of the  page at hyperphysics that you link to have left out an important step; they have not considered exactly how that highly oscillatory probability (the envelope of which looks classical) gets smeared out to remove the oscillations...
And, I'll bet that now you can guess the answer as to how that happens...
