Liouville equation solution $\rho (q,p,t)$ for 1D-oscillator

I'd like to consider the 1D-oscillator phase space probability density evolution problem with the ordinary Hamiltonian $$H = \frac{p^2}{2m} + \frac{kq^2}{2}.$$

Then the Liouville theorem is just

$$\frac{{\partial \rho }}{{\partial t}} = { {\frac{{\partial \rho }}{{\partial {q}}}\frac{{\partial H}}{{\partial {p}}} - \frac{{\partial \rho }}{{\partial {p}}}\frac{{\partial H}}{{\partial {q}}}}}$$

for the case of one particle.

Assume that the initial distribution of $$q$$ and $$p$$ for a particle is normal
$$\rho (q,p,t = 0) = \frac{1}{{2\pi }}{e^{ - {q^2} - {p^2}/2}}$$

Now we can find the Hamiltonian partial derivatives $$\frac{\partial H}{\partial p} = p/m$$ and $$\frac{\partial H}{\partial q} = -kq$$.

After this step I'm confused how we can find the evolution of the probability?

How are $$p$$, $$q$$ related with distribution $$\rho (q,p,t = 0)$$?

What is the $$\frac{{\partial \rho }}{\partial {q}}$$ and $$\frac{{\partial \rho }}{\partial {p}}$$ when we even don't know $$\rho (q,p,t)$$ in arbitrary moment of time $$t$$? Maybe my undestanding is wrong, but I just want to know $$\rho (q,p,t)$$.

Consider a set of $$N$$ harmonic oscillators. $$\rho(q,p,t=0)$$ gives the probability to find one of these oscillators in the state $$(q,p)$$ at time $$t$$. Using the equations of motion (Hamilton equations in your case) $${dq\over dt}={p\over m},\hskip 1cm {dp\over dt}=-kq$$ you can tell that if there was an oscillator in the state $$(q_0,p_0)$$ then this oscillator is in the state \left\{\eqalign{ q(t)&=q_0\cos\omega t+{p_0\over\omega}\sin\omega t\cr p(t)&=-\omega x_0\sin\omega t+p_0\cos\omega t\cr }\right. at time $$t$$. Equivalently, you can tell that, if there is an oscillator in the state $$(q,p)$$ at time $$t$$, then it was in the state \left\{\eqalign{ q_0&=q\cos\omega t-{p\over\omega}\sin\omega t\cr p_0&=\omega q\sin\omega t+p\cos\omega t\cr }\right. at time 0. Therefore, the probability to find an oscillator in the state $$(q,p)$$ at time $$t$$ is $$\rho(q(t),p(t),t)=\rho(q(0),p(0),0)$$ which reads $$\rho(q,p,t)=\rho\Big(q_0=q\cos\omega t-{p\over\omega}\sin\omega t, p_0=\omega q\sin\omega t+p\cos\omega t,0\Big)$$ This is actually the solution of Liouville equation $${d\rho\over dt}={\partial\rho\over\partial t} +\dot q{\partial\rho\over\partial q} +\dot p{\partial\rho\over\partial p}=0$$ and this way of solving it is known as the method of characteristics.

• Do I just need to put the expressions for $q(t)$ and $p(t)$ into the initial Gauss distribution $\rho (q,p,t = 0) = \frac{1}{{2\pi }}{e^{ - {q^2} - {p^2}/2}}$? Great! Commented May 19, 2020 at 13:51
• Yes : $\rho(q,p,t)={1\over 2\pi}e^{-[q\cos\omega t-{p\over\omega}\sin\omega]^2 t-{1\over 2}[\omega q\sin\omega t+p\cos\omega t]^2}$. Commented May 19, 2020 at 13:56
• For dimensional and physical reasone, I guess that our initial distribution is ${1\over 2\pi}e^{-{1\over 2}[p^2+\omega^2q^2]}$, no ? Then, since energy is conserved by the dynamics, you get $e^{-H(0)}=e^{-H(t)}$. Commented May 19, 2020 at 14:00
• Yeah, I think you're correct. Thank you Commented May 19, 2020 at 14:08

According to Liouville's theorem,

$$\frac{d\rho}{dt}=\frac{\partial \rho}{\partial t}+\{\rho,H\}=0$$ where {,} stand for the Poisson brackets. Now, in equilibrium statistical mechanics, we are concerned about the behaviour of systems in equilibrium when $$\rho$$ has no explicit time dependence, i.e., when$$\frac{\partial \rho}{\partial t}=0 .$$ So, to make these two equations valid simultaneously (or in other words, for the Liouville's theorem to hold when the system is in equilibrium), we need to have $$\{\rho,H\}=0$$ . This is possible if $$\rho(q,p)$$ is equal to some constant or more generally if $$\rho=\rho[H(q,p)]$$ , i.e., if $$\rho$$ has its dependence on q and p only through $$H$$.

To see this, you may read Statistical Mechanics by R. K. Pathria and Paule D. Beale (Second Chapter of Third Edition).

Now, the form of $$\rho$$ depends on the type of the ensemble concerned; for example, $$\rho$$=constant for a microcanonical ensemble because all the microstates are equally probable here. In a canonical ensemble, $$\rho$$ is proportional to $$e^{-H(X)/k_BT}$$ where $$X=\{q_i,p_i\}$$; $$i=\{1,2,...,N\}$$ , for N particles in one dimension.

For your problem, I do not see how the initial distribution that you have mentioned is relevant. However, if we assume that there is an isolated system of energy U and volume V comprising of N non-interacting 1D linear harmonic oscillators, then it forms a microcanonical ensemble where the suitable form for $$\rho(q,p)$$ is as follows :
$$\rho(q,p)=\frac{\delta(H_N(X)-U)}{\Sigma}$$ where $$\Sigma=\displaystyle{\int}_{\Gamma}\frac{d^{2N}X}{h^N}\delta(H_N(X)-U)$$ is the microcanonical partition function. Here, $$H_N(X)=\displaystyle{\sum_{i=1}^N}\frac{p_i^2}{2m}+\frac{kq_i^2}{2}$$ . So,
$$\Sigma=\displaystyle{\int}_{\Gamma}\frac{d^{2N}X}{h^N}\delta(H_N(X)-U)= \displaystyle{\int_{U{\le}H_N(X){\le}U+\Delta}}\frac{d^{2N}X}{h^N}$$ $$=\displaystyle{\int_{U{\le}\sum_{i=1}^N\frac{p_i^2}{2m}+\frac{kq_i^2}{2}{\le}U+\Delta}}\frac{d^{2N}X}{h^N} ={(2m)}^{N/2}{\left(\frac{2}{k}\right)}^{N/2}\displaystyle{\int_{U{\le}\sum_{i=1}^N(p_i^2+q_i^2){\le}U+\Delta}}\frac{{(dpdq)}^N}{h^N}$$
$$=\frac{{(2m)}^{N/2}}{h^N}{\left(\frac{2}{k}\right)}^{N/2}\times\,Volume\,of\, a \,shell\, between\, the\, 2N-dimensional\, hyperspheres\, of\, radii\, \sqrt{U}\, and\, \sqrt{U+\Delta}\,\,[where\, \Delta<
$$=\frac{{(2m)}^{N/2}}{h^N}{\left(\frac{2}{k}\right)}^{N/2}\times\,surface\,of\,a\,2N-dimensional\,sphere\,of\,radius\, \sqrt{U}\,\times\,thickness\,of\,the\,shell$$
$$=\frac{{(2m)}^{N/2}}{h^N}{\left(\frac{2}{k}\right)}^{N/2}\frac{(2N){\pi}^{N}({\sqrt{U}}^{2N-1})}{\Gamma(N+1)}\times\frac{\Delta}{2\sqrt{U}}$$
which you can further simplify.
This is for distinguishable oscillators. If they are indistinguishable, $$\Sigma$$ should be divided by $$N!$$ . Putting the expression for $$\Sigma$$ in that of $$\rho$$, you will get the desired expression.

• It's quite helpful, thanks Commented May 19, 2020 at 15:14