Does quantum mechanics ever really quantize classical mechanics?  I was curious about a physics question which I thought might be suitable for mathoverflow. I looked at the answer to this question, but it's not what I'm looking for.
Basically, classical mechanics and the $\hbar \to 0$ limit of quantum mechanics study the action of the same algebra on very different representations. I'm curious whether there is a good physical explanation for why as you degenerate to the $\hbar \to 0$ limit the algebra of observables degenerates to the same Poisson algebra appearing in classical mechanics, but the relevant representation changes significantly. Specifically, (non-relativistically) classical systems have evolution $\frac{d}{dt} \rho = \{H, \rho\}$ and quantum systems have evolution $\frac{d}{dt}\psi = [H, \psi]$ (up to some constants). So far these look analogous, but in classical mechanics, the density function $\rho$ is itself a function on the phase space (i.e. a vector of the regular representation), whereas in quantum mechanics, $\psi$ is a is just a function (or something) on the $x_i$ themselves - i.e. a vector of a representation of "square-root dimension" (half-dimensional singular support)! 
My guess is that this is a many-particle phenomenon, and a fully honest answer to "why do we observe classical mechanics" will probably involve a serious study of deconherence and questions of "what is observation", etc. 
But I'm curious if there is a heuristic way to see why the algebra that's acting is the same (and in what way the representation is allowed to change: e.g., is there some embedding of the regular representation in a tensor product of irreducible ones?)
 A: It is perhaps helpful to distinguish between four types of mechanics here:


*

*Pure-state classical mechanics.  Here, the mechanics are classical, and the system is described by a single point $(q,p)$ in phase space.  This point evolves via Hamilton's equations of motion $\partial_t q = \frac{\partial H}{\partial p}; \partial_t p = - \frac{\partial H}{\partial q}$.

*Mixed-state classical mechanics.  Here, the mechanics are classical, and the system is described by a probability density function $\rho(q,p)$ on phase space (this density may be a generalised function, e.g. a Dirac delta, rather than a classical function).  This density function evolves via the advection equation $\partial_t \rho = \partial_p ( \rho \partial_q H ) - \partial_q (\rho \partial_p H ) = \{H,\rho\}$.

*Pure-state quantum mechanics.  Here, the mechanics are quantum, and the system is described by a wave function $|\psi\rangle$ in a Hilbert space.  This wave function evolves via Schrödinger's equation of motion $\partial_t |\psi \rangle = \frac{1}{i\hbar} H |\psi\rangle$.

*Mixed-state quantum mechanics.  Here, the mechanics are quantum, and the system is described by a density matrix $\rho$ (a positive semi-definite trace one operator on a Hilbert space).  This density matrix evolves by the von Neumann evolution equation $\partial_t \rho = \frac{1}{i\hbar} [H,\rho]$.


In both the classical and quantum regimes, a mixed state can be viewed as a convex (or classical) superposition of pure states (with a pure classical state $(q,p)$ identified with the Dirac probability density function $\delta_{(q,p)}$, and a pure quantum state $|\psi \rangle$ identified with a pure density matrix $|\psi \rangle \langle \psi|$).  So in principle the pure-state mechanics describes the mixed-state mechanics completely (albeit with the caveat that in the quantum case, in contrast to the classical case, the decomposition of a mixed state as a superposition of pure states is non-unique).  However, the correspondence principle is clearest to see at the mixed state level, i.e. to compare 2. with 4. in the semiclassical limit $\hbar \to 0$, rather than comparing 1. with 3..  Indeed, any density matrix $\rho$ has a Wigner transform $\tilde \rho$, which is a function on phase space defined via duality as $\int \tilde \rho(q,p) A(q,p)\ dq dp = \hbox{tr}( \rho \hbox{Op}(A) )$ for any classical observable $A$, where $\hbox{Op}(A)$ is the (Weyl) quantisation of $A$ (i.e. the Wigner transform is the adjoint of the quantisation operator).  This Wigner transform $\tilde \rho$ will usually not be non-negative, and hence will not be a classical probability density function, but in semiclassical regimes it is often the case that $\tilde \rho$ will tend (in a suitable weak sense) to a classical probability density when $\hbar \to 0$, which will then evolve by the classical advection equation.  This is the dual to the assertion that the quantum Heisenberg equation $\partial_t A = \frac{i}{\hbar} [H,A]$ for the evolution of quantum observables converges to the classical Poisson equation $\partial_t A = -\{ H,A\}$ for the evolution of classical observables in the semiclassical limit $\hbar \to 0$.
There is still a correspondence at the level of 1. and 3., but it is a bit trickier to see; one has to restrict to things like "gaussian beam" type solutions $|\psi \rangle$ to the Schrödinger equation that are well localised in both position and momentum space, in order to get a classical limit that is a pure state rather than a mixed state.  (An arbitrary wave function would instead get a "phase space portrait" which in the semiclassical limit becomes [assuming some equicontinuity and tightness, and possibly after passing to a subsequence, as noted in comments] a mixed state from 2., rather than a pure state from 1.).
 
A: I interpret your question as a query into a mathematical formulation of quantum decoherence, which is the process by which a partial trace of the quantum mechanical density operator $\hat\rho$ reduces to the classical phase space function $\rho$.
A simple case where this process can be analyzed in much mathematical detail is described in Decoherence in a Two-Particle Model (2001): We consider a simple one dimensional quantum system consisting of a heavy and a light particle interacting via a point interaction. The initial state is chosen to be a product state, with the heavy particle described by a coherent superposition of two spatially separated wave packets with opposite momentum and the light particle localized in the region between the two wave packets. We characterize the asymptotic dynamics of the system in the limit of small mass ratio, with an explicit control of the error. We derive the corresponding reduced density matrix for the heavy particle and explicitly compute the (partial) decoherence effect for the heavy particle induced by the presence of the light one for a particular set up of the parameters.
For the algebraic aspects, see
An Algebraic Formulation of Quantum Decoherence: An algebraic formalism for quantum decoherence in systems with continuous evolution spectrum is introduced. A certain subalgebra, dense in the characteristic algebra of the system, is defined in such a way that Riemann-Lebesgue theorem can be used to explain decoherence in a well defined final pointer basis.
A: There exists an interesting operatorial (Hilbert space) approach to classical mechanics pionereed by Koopman and von Neumann: http://arxiv.org/abs/quant-ph/0301172 (Topics in Koopman-von Neumann Theory, by D. Mauro). 
Within this formalism (or, more precisely, in its  functional integral version), quantization is mysteriously associated to the freezing to zero of two Grassmannian partners of time: http://arxiv.org/abs/quant-ph/0308101 (Time and Geometric Quantization, by A.A. Abrikosov Jr, E. Gozzi and D. Mauro).
A: I thought it might be nice to couple Terry Tao's great general answer by showing we can write down an explicit limit to the classical case for the simple harmonic oscillator. These solutions are an example of "coherent states". I learned this from an old blog post by John Baez which I can't find right now; Wikipedia has a less helpful exposition.
We work with an oscillator of frequency $\omega$, so the potential energy is $(1/2) m \omega^2 x^2$ and Schrödinger's equation is
$$i \hbar \frac{\partial}{\partial t} \psi = - \frac{\hbar^2}{2 m} \frac{\partial^2}{(\partial x)^2} \psi + \frac{m \omega^2}{2} x^2 \psi.$$
As usual, it is convenient to define
$$a = \sqrt{\frac{m \omega}{2 \hbar}}\left(x+\frac{\hbar}{m \omega} \frac{\partial}{\partial x} \right) \quad \mbox{annihilation}$$
$$a^{\dagger} = \sqrt{\frac{m \omega}{2 \hbar}}\left(x-\frac{\hbar}{m \omega} \frac{\partial}{\partial x} \right) \quad \mbox{creation}.$$
The lowest energy state is the kernel of $a$, namely $\psi_0 := \exp(-m \omega x^2/(2 \hbar))$; it gives rise to the solution $e^{i \omega t/2} \psi_0$. Then $\frac{1}{\sqrt{n!}} (a^{\dagger})^n \psi_0$ is the $n$-th energy state, so $e^{i (n+1/2) \omega t} (a^{\dagger})^n \psi_0$ is the $n$-th solution to the time dependent equation (up to normalization). I prefer to rewrite this as $(e^{i \omega t} a^{\dagger})^n (e^{i \omega t/2} \psi_0)$.
If $F(z)=\sum f_n z^n$ is any power series then, at least formally, $F(e^{i \omega t} a^{\dagger})(e^{i \omega t/2} \psi_0)$ is a solution of Schrödinger's equation, since it is a linear combination of the pure energy states above. 
In particular, take $F(z) = \exp(C z)$ for some scalar $C$. So
$$\exp\left( C e^{i \omega t} \sqrt{\frac{m \omega}{2 \hbar}}\left(x-\frac{\hbar}{m \omega} \frac{\partial}{\partial x} \right) \right) (e^{i \omega t/2} \psi_0)$$
solve Schrödinger's equation. We reduce clutter by setting $C\sqrt{\frac{\hbar}{2 m \omega}}=R$; the constant $R$ has units of distance.
So our solution is 
$$\exp\left( e^{i \omega t} \left(\frac{R m \omega}{\hbar} x-R \frac{\partial}{\partial x} \right) \right) (e^{i \omega t/2} \psi_0)$$
Now, the commutator of $\tfrac{R m \omega}{\hbar} x$ and $R \tfrac{\partial}{\partial x}$ is $\tfrac{m R^2 \omega}{\hbar}$, which commutes with both $x$ and $\partial/\partial x$. So, by Baker-Cambell-Hausdorff (and discarding some global constants) we can rewrite $(\ast)$ as
$$\exp(e^{2 i\omega t} \tfrac{m R^2 \omega}{\hbar}) \exp(e^{i \omega t}\tfrac{R m \omega}{\hbar} x ) \exp\left( -e^{i \omega t} R \frac{\partial}{\partial x} \right) (e^{i \omega t/2} \psi_0).$$
The exponential of differentiation is translation, so this is
$$\exp(e^{2 i\omega t} \tfrac{m R^2 \omega}{\hbar}) \exp(e^{i \omega t}\tfrac{R m \omega}{\hbar} x )  (e^{i \omega t/2} \psi_0(x-R e^{i \omega t})).$$
One can then do a bunch of work translating each formula into its real and complex part, which I omit. At the end of the day, one gets a solution to Schrödinger's equation which roughly looks like
$$e^{i A(t)} \exp\left(\tfrac{m \omega}{\hbar} \left[-(x-R \cos(\omega t))^2/2 - i R \sin(\omega t) x \right] \right).$$
Here $A$ is a big messy function I am unwilling to work out.
This is the sort of gaussian beam solution Terry was talking about -- it is localized both in position and in Fourier space. In position space, it is a Gaussian centered at $x=R \cos (\omega t)$. As $\hbar \to 0$ (with $R$ fixed), the Gaussian becomes tighter and tighter until, in the limit, it is a delta function at $R \cos (\omega t)$ -- the classical solution to the problem. Meanwhile, the momentum is a Gaussian centered at $-m R \omega \sin(\omega t)$. Again, as $\hbar \to 0$, the Gaussian becomes a delta function at $-m R \omega \sin(\omega t)$ -- the classical solution. (Of course, you could ignore all the discussion about $a$ and $a^{\dagger}$ and just directly check that this solves Schrödinger's equation. If you do, please let me know what constants I left out!)
If one tries to take the $\hbar \to 0$ limit of some simpler solutions like the pure energy states, they bunch up at the origin while spreading out over all of momentum space. You need a moderately complicated solution like this to get both limits to make sense.
I will close by noting a heuristic way to think about $\exp(C a^{\dagger})$. The coefficient of the $n$-th energy state is $C^n/\sqrt{n!}$ (putting in the correct normalization constant.) So, if we observe the energy of this particle, we have probability proportional to $C^{2n}/n!$ of getting the answer $(n+1/2) \hbar \omega$. In other words, the energy of this particle is a Poisson random variable with expected value $C^2 \hbar \omega+\hbar \omega/2$. Plugging in for $C^2$, this is $m R^2 \omega^2/2+\hbar \omega/2$. The $m R^2 \omega^2/2$ term is the energy of the classical solution. So this solution may be thought of as the best attempt to mimic an energy of $m R^2 \omega^2/2$ when we only have access to the discrete levels $n \hbar \omega$.
A: It might be worth to add, that Planck found "his" constant $h$ by considering a many-body-problem indeed. 
At this time, it was quite a standard method in statistical mechanics, to divide the physical phase space in discrete small phase volumes ("cells") and to count them. So, statistical mechanics at this time already contained a double limit procedure, by letting the number of particles $N \to \infty$ (and the total volume $V \to \infty$ as well, to obtain a thermodynamic limit), and by letting the discretizing phase space volume unit $\nu \to 0+$. 
From a today's point of view, to obtain classical many-body statistical mechanics, one first takes $N \to \infty$ in the thermodynamic limit of quantum statistical mechanics, and then $\hbar \to 0+$ afterwards to obtain the (semi-)classical limit.
This is in a sense, and mathematically spoken from a today's point of view, what Planck found out (without saying so).
Classical statistical mechanics would have been obtained by first taking the "(semi-)classical limit" $\hbar \to 0+$, and the thermodynamics limit $N \to \infty$ afterwards.
A nice reference is Longair: Theoretical Concepts in Physics, 2ndEd, 2003, chap 13 (see also chap 12).
One also has to understand, that physics is not a formal science, like mathematics, but its ultimate goal is to understand real phenomena. There's no one-particle-system anywhere in the world. There exist only many-body-systems. Considering a system with one particle only or with one degree of freedom etc. stems from big reduction procedures and perturbation theory.
John von Neumann's Quantum Mechanics' book has a chapter about the coupling of the observed quantity to the environment. 
In elemetary form, the procedure to obtain "physically treatable" models is as follows. One considers a certain part of reality an tries to specify it by declaring a boundary, separating the "observed system" and the "environment". One takes a thermodynamic limit for the environment. Then we can speak about intensive quantities such as temperature. These quantities determine a (one or maybe several) equilibrium state of the observed system. We make a perturbation ansatz for the observed system by considering states near the equilibrium. We try to get rid of all detailed interactions of the observed system with the detailed "particles" of the environment, by considering only averaged interactions with the environment. In general, the obtained equations for a subsystem would be non-local (both in space and time). We localize near the equilibrium to obtain equations for our observed subsystem. The data of the environment enter only through some (not too many) parameters. 
Here is where most models of physics start.
These considerations make also clear, why one has to take into account mixed states (both in quantum statistical mechanics and in classical statistical mechanics) and not only pure states. The many averaging procedures mentioned above will in general provide a mixed state for the observed subsystem, not a pure state.
Pure states are mainly for helping in modelling, and for analyzing some phenomena - because mixed states can in general be written as a combination of pure states.
I would not consider quantum mechanics to "quantize" classical mechanics, because there is no unique procedure classical -> quantum. It's rather the other way round, that quantum mechanics models have also a (semi-)classical limit.
It would be physically misleading to think, that geometric quantization or deformation quantization are the answer to the question of how to model a classical system quantum mechanically. That's basically not true. The quantization procedures help only in identifying which kind of quantum models could be appropriate, and which are clearly to be discarded. 
The real procedure in physics is, to guess a quantum mechanical model (most often by experience, and as a combination of simpler models that are already understood) and to look by experiment how good it describes reality. Classical mechanics only helps in guiding modelling eg. via geometry (for instance, a geometry of indistinguishable fermions is not so easy to visualise...).
Physicists' thinking is in some sense opposite to the thinking of mathematicians. A mathematician usually has a system of axioms for the "general case" and the specializes by making further assumptions.
In physical modelling for new situations, there are basically no axioms (except very general ones), but there are very specific models, that have been tested in specific situations. These models are then combined, like in a modular construction system. That's also, why physics textbooks almost never have axioms, but a plethora of examples and exercises to be calculated. Physical modelling in this sense is "bottom up", where the method of how a mathematician models is rather "top down". Mathematicians model only in situations, where the physics has been settled, eg. take the already know general equations of continuum mechanics and specifiy by choosing appropriate material constants (or functionals). Physical modelling usually occurrs in situations, where the "general equations" are unknown, so they can't be specialized to the specific situation.
