mathematical physics - Probability density in Hamiltonian Mechanics

I am currently studying Liouville's theorem compare wikipedia

and there this mysterious probability density $\rho$ appears and I was wondering how one can determine this quantity analytically for a given problem.

Let's think of the harmonic oscillator as a primitve example, how can I get an analytical representation of $\rho$?

Answer

Liouville's theorem describes the time-evolution of the density, in canonical ensemble where the energy thereof is conservative, of system points in phase space. A common example is a closed system with Hamiltonian $\mathcal{H}=\frac{1}{2}p^2+V(q)$. According to Liouville's equation, the probability density distribution $\rho(q,p;t)$ of phase space remains invariant along trajectory. That is $$\frac{d\rho}{dt}=\frac{\partial\rho}{\partial t}+\frac{\partial\rho}{\partial q}\frac{\partial\mathcal{H}}{\partial p}-\frac{\partial\rho}{\partial p}\frac{\partial\mathcal{H}}{\partial q}=\frac{\partial\rho}{\partial t}+\frac{\partial\rho}{\partial q}p-\frac{\partial\rho}{\partial p}\nabla V=0$$ It's a linear first-order pde and we can solve it by method of characteristics. The characteristic equation is $$dt=\frac{dq}{p}=-\frac{dp}{\nabla V}$$ which gives two first integrals. $$\mathcal{H}=\frac{1}{2}p^2+V(q)$$ $$q_0=q-\int pdt$$ Thus the solution is $$\rho(q,p;t)=F(\mathcal{H},q-\int pdt)$$ for any function bivariate $F(x,y)$.

Assume the initial distribution is multivariate normal distribution $\mathcal{N}({\bf0,\mathcal{I}})$, namely $$\rho(q,p;0)=F(\mathcal{H},q_0)=(2\pi)^{-\frac{n}{2}}\exp\{-\frac{1}{2}q_0^Tq_0\}$$ Then the final solution $$\rho(q,p;t)\propto\exp\{-\mathcal{H(q(t),p(t))}\}$$