Does quantum mechanics ever really quantize classical mechanics?

It is perhaps helpful to distinguish between four types of mechanics here:

1. 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}$.
2. 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\}$.
3. 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$.
4. 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.).

 

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$.

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.