# 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$?

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))}\}$$

• Plug in V(q) being the attractive coulomb potential. If you then integrate(just let me say in x-direction) $\int_{-\infty}^{\infty} \rho(x,y,z,p,t) dx$ this integral diverges. That does not seem reasonable to me. – Xin Wang Sep 3 '13 at 10:59
• There's no problem since you have to guarantee Hamiltonian is finite. That is, if you assume Coulomb potential, then initial distribution cannot be normal – Shuchang Sep 3 '13 at 11:55
• mhmm, new question: does the concept of a probability density distribution apply to systems with only one particle? I mean, sure, there are applications of this concept in statistical mechanics, but if I want to solve a mechanical problem, where I am dealing with two charges(let's stay with the coulomb problem) and both initial values(twice the position and velocities are known). does the concept of a probability density apply to this case? In that case, my probability density would be rather something like a delta function at all states, right? – Xin Wang Sep 3 '13 at 11:59
• Yes it is. But if you handle conservative systems of only few particles, there's no need to use Liouville's theorem. Hamiltonian mechanics works well. – Shuchang Sep 3 '13 at 12:03
• ah, then i got your idea – Xin Wang Sep 3 '13 at 12:28

The probability density is the state of the system in this setting. While the hamiltonian contains all the information about how the system is set up, what forces are acting, and so on, the probability density contains all the physically measurable information about whereabouts in phase space your system actually is.

This is a complete shift in paradigm to other approaches to mechanics. In quantum mechanics, for example, you have a fundamental state - the wavefunction, $\psi$ - and you find $\rho=|\psi|^2$ in terms of it. In hamiltonian mechanics you might play the game of forgetting when you started the clock, and seeing where you might find the particle most often (see, e.g. "Quantum and classical probability distributions for position and momentum," RW Robinett, Am. J. Phys. 63 no. 9, p. 823 (1995), available here). This is different. Since you don't know where you've put your system, or where it might be now, the probability density $\rho$ is your fundamental handle on the state of the system, and you can't trade it for other objects.

Thus, given a specific initial state (and a specific hamiltonian), you have to solve Liouville's equation to obtain $\rho$ for later times. Depending on the system and the initial state this may or may not be possible analytically, or you may even find yourself in numerical difficulties. In general, though, this is relatively hard, since your equation of motion has changed from a system of ODEs (as in hamiltonian mechanics) to a relatively complex partial differential equation.

That said, the harmonic oscillator is in fact a very, very special case for liouvillian mechanics. This is because of its harmonicity: more precisely, because all the orbits have the same period. This implies that phase space moves rigidly, rotating around the equilibrium position at the oscillation frequency, and it carries the probability density with it.

Thus, for a harmonic hamiltonian $H=\frac12(p^2+q^2)$, the solution to Liouville's equation $$\left[\frac\partial{\partial t}+p\frac\partial{\partial q}-q\frac\partial{\partial p} \right]\rho=0$$ under an arbitrary initial condition $\rho_0(q,p)$ is easily found to be $$\rho(q,p,t)=\rho_0(q\cos(t)-p\sin(t),q\sin(t)+p\cos(t)).$$ (You should, of course, do the corresponding calculation! For a good exercise, add units and check that everything works out.) It should be pretty clear that what this solution does is trace the classical trajectory back from the point $(q,p)$ at time $t$ to where it would have started in at time 0, and retrieve the original particle density there.

I imagine there is an equivalent method to obtain (at least formal) solutions to Liouville's equations in terms of the corresponding hamiltonian trajectories. (A brief foray into the web didn't turn it up, but I'm pretty sure it exists.) Thus you'd reduce the hard PDE problem to the ODEs it was originally constructed from. Do note, though, that in general things will be at least slightly more complicated than the above, because if the phase-space motion isn't rigid - if trajectories bunch up or spread out - the probability density will correspondingly increase or decrease.

• @Lipschitz I don't understand either comment. The probability density is indeed a scalar; I don't see where you infer otherwise. The integral will always converge if it converges at $t=0$, because Liouville's equation conserves probability (i.e. it can be written as $\frac\partial{\partial t}\rho+\nabla\cdot\vec j=0$). – Emilio Pisanty Sep 3 '13 at 10:53
• well you have a comma in your equation, so this does not look like a scalar. – Xin Wang Sep 3 '13 at 11:00
• The comma is there as in all functions of two variables, like $f(x,y)$, to separate the different variables. $\rho_0=\rho_0(q_0,p_0)$ is a function of $q_0$ and $p_0$, and here $q_0=q\cos(t)-p\sin(t)$ and $p_0=q\sin(t)+p\cos(t)$. Hence the comma. – Emilio Pisanty Sep 3 '13 at 11:02
• ah, sorry this was unclear to me since you started with $\rho_0(q,p)$ first rather than $\rho_0(q_0,p_0)$. In that case I have just one question lest: how did you determine $q_0 = q cos(t)- psin(t)$ and $p_0 = q sin(t)+pcos(t)$. I guess this comes somehow from the solution of the mechanical problem, right? In this case, it should be always possible to do it that way. – Xin Wang Sep 3 '13 at 11:11

Consider some classical mechanical system with canonical coordinates and momenta $(p,q)$ on its phase space $\Gamma$. We can imagine building a large number of identical copies of this system. We call such a set of copies an ensemble.

Next, we can prepare the initial conditions of each element of the ensemble however we like. If we imagine that the state of each element of the ensemble is represented as a point in $\Gamma$, then the whole ensemble will look like a cloud of points in phase space. If we imagine that the number of points is very large, then we can imagine representing the cloud in terms of a density on phase space. If, for example, all of the members of the ensemble are initially prepared in nearly the same state, then the cloud will look like a clump around the corresponding point in phase space.

Once we specify the initial behavior of the ensemble however, the Hamiltonian evolution will determine how the cloud evolves with time. In fact, the equation that governs the evolution of the phase density is \begin{align} \frac{\partial\rho}{\partial t} = \{H,\rho\} \end{align}

Constructing an example is easy. Simply choose your favorite probability density on phase space at the initial time, and use the evolution equations to obtain it for later times!

• sorry, but after all what i have seen, is the idea to determine the density probably the hardest approach to solve a mechanical problem. so i am just asking for one simple system to be solved.(if you do not like the harmonic oscillator, you may pick something else) – Xin Wang Sep 3 '13 at 8:18
• @Lipschitz The explicit construction for the simple harmonic oscillator is given in Emilio's answer. As both he and I remark, you have the freedom to choose $\rho_0$, the initial phase density. Such a choice corresponds to how the statistical state of the system is prepared initially. – joshphysics Sep 3 '13 at 8:26
• and I argued why I doubt that this is a solution.(Shuchang Zhang ended up with a completely different solution) – Xin Wang Sep 3 '13 at 10:53
• @Lipschitz That is because Shuchang chose the initial distribution to be a normal distribution. As indicated by both Emilio's answer, and my answer, (and implicitly in Shuchang's answer) you can choose the initial distribution, so there is an infinite number of possible answers, each of which is correct. – joshphysics Sep 3 '13 at 15:12