Liouville's Theorem

Definition

Picture not one system but a whole ensemble of identical systems, each launched from a slightly different starting condition. In phase space — the abstract space whose axes are all the positions q and all the momenta p — that ensemble is a cloud of points with some density ρ(q, p, t).

Liouville's theorem (Joseph Liouville, 1838) states that ρ is conserved along the flow:

dρ/dt = 0     (the total / convective derivative)

Ride along with any single system as it moves through phase space, and the density of neighbours around you never changes. Geometrically, this means the volume occupied by a chunk of the cloud is conserved as the chunk is carried along by the dynamics. The blob can shear, stretch into a long thin ribbon, and fold over on itself like taffy — but its total area (in 2 dimensions) or 2N-dimensional volume never grows or shrinks. Phase-space volume is incompressible.

How it works

The result follows from two facts glued together. First, ensemble members are neither created nor destroyed, so ρ obeys a continuity equation in phase space:

∂ρ/∂t + ∇·(ρ v) = 0     where v = (q̇, ṗ) is the phase-space velocity

Second, the flow is generated by Hamilton's equations:

q̇ = ∂H/∂p        ṗ = −∂H/∂q

Compute the divergence of v. For each degree of freedom i:

∂q̇ᵢ/∂qᵢ + ∂ṗᵢ/∂pᵢ
  = ∂/∂qᵢ (∂H/∂pᵢ) + ∂/∂pᵢ (−∂H/∂qᵢ)
  = ∂²H/∂qᵢ∂pᵢ − ∂²H/∂pᵢ∂qᵢ
  = 0      (mixed partials commute)

So ∇·v = 0: the Hamiltonian flow is divergence-free, exactly like an incompressible fluid. Expanding the continuity equation with ∇·(ρv) = ρ(∇·v) + v·∇ρ and dropping the zero term:

∂ρ/∂t + v·∇ρ = 0     ⟹     dρ/dt = ∂ρ/∂t + {ρ, H} = 0

where {ρ, H} is the Poisson bracket. The convective derivative of ρ vanishes — that is Liouville's theorem. Because the flow is volume-preserving, the Jacobian of the time-evolution map from t to t+Δt has determinant exactly 1; this is the statement that phase volume is conserved.

A worked example with numbers

Take the simplest non-trivial Hamiltonian system, a particle moving freely along a line: H = p²/2m, with no force. Hamilton's equations give q̇ = p/m and ṗ = 0. A particle's momentum is frozen; its position drifts at constant speed p/m.

Now seed an ensemble as a rectangle in phase space at t = 0: positions spread over Δq₀ = 2 m, momenta over Δp₀ = 4 kg·m/s. Initial area:

A₀ = Δq₀ · Δp₀ = 2 × 4 = 8 m·(kg·m/s)

Each point evolves as q(t) = q₀ + (p₀/m)·t. Fast points (high p) outrun slow ones, so the rectangle shears into a parallelogram tilted along q. After t = 3 s with m = 1 kg, the top edge has slid forward by (4/2)/1 × 3 — the spread in q grows because the momentum spread is now expressed as a position spread. The horizontal extent stretches, yet:

Area of a sheared parallelogram = base × height = Δq₀ × Δp₀ = 8

The shear transformation has determinant 1, so the area is still 8 — unchanged. The blob got longer and thinner in one diagonal direction, but the bookkeeping always balances: every metre the cloud spreads in position, it must have been "paid for" by an existing spread in momentum. That is incompressibility in action, with concrete units. Run the same ensemble through a harmonic oscillator instead and the rectangle rotates rigidly in phase space; through a chaotic system and it folds into a fractal filament — but in all three the measured area stays 8.

Foundation of the microcanonical ensemble

Why does any of this matter? Because it is what makes statistical mechanics consistent. The microcanonical ensemble — the bedrock of the whole subject — assumes that for an isolated system at fixed energy E, every accessible microstate on the energy shell is equally likely. In density language, ρ is uniform on the shell H = E and zero elsewhere.

For that assumption to mean anything, the uniform density must stay uniform as time passes. If phase volume were not conserved, an initially flat ρ would pile up in some regions and drain out of others within one dynamical time, and "equal a-priori probability" would be a fiction by the next instant. Liouville's theorem rescues it: a density that is a function only of conserved quantities, ρ = ρ(H), satisfies {ρ, H} = 0, so

∂ρ/∂t = −{ρ, H} = 0     ⟹     ρ is stationary

The uniform-on-the-shell distribution is therefore an equilibrium distribution — it is invariant under the exact dynamics. The same argument, with ρ ∝ e^(−H/kT), produces the stationary canonical ensemble. Liouville's theorem is the silent assumption underneath every partition function you have ever written.

Variants and regimes

System / framework	Conserved object	Volume behaviour	Governing equation
Hamiltonian (frictionless)	Fine-grained ρ	Exactly preserved	dρ/dt = 0
Symplectic integrator (numerics)	Phase volume (machine-precision)	Preserved by construction	Discrete symplectic map, det = 1
Dissipative (drag −γp)	None — attractors	Contracts as e^(−γt)	∇·v = −γ < 0
Stochastic (noise + drag)	Probability (not volume)	Spreads / settles to steady state	Fokker–Planck equation
Quantum (closed)	Eigenvalues of ρ̂, entropy	Preserved by unitarity	von Neumann: iℏ∂ρ̂/∂t = [Ĥ, ρ̂]
Coarse-grained ensemble	Total probability only	Effective volume grows (mixing)	Boltzmann / master equation
Relativistic phase space	Invariant d³x d³p	Lorentz-invariant, preserved	Relativistic Vlasov / Boltzmann

Common pitfalls and misconceptions

• "Conserved volume means the cloud keeps its shape." No. The shape can deform without bound — into ribbons, spirals, fractal whorls. Only the measure is fixed. The visual drama of Liouville flow is precisely that the blob distorts violently while its area sits perfectly still.
• "Liouville's theorem forbids entropy increase." It forbids fine-grained entropy from changing. Real entropy growth comes from coarse-graining: the filaments get thinner than any measurement scale, so the smeared-out density fills the shell and coarse-grained entropy rises. The two statements coexist.
• "It applies to any equation of motion." Only to Hamiltonian (or symplectic, divergence-free) flows. Add friction, a thermostat, or a non-conservative force and volume is no longer preserved. This is a feature, not a bug — it is why dissipative systems can have attractors.
• "Volume conservation lets me focus a beam arbitrarily." The opposite. Squeeze the position spread Δq and the momentum spread Δp must grow to keep ΔqΔp fixed. Accelerator physicists call this conserved quantity the emittance; you can rotate it between planes but never beat it with linear optics. This is the classical shadow of the uncertainty principle.
• "ρ is the probability of one system." ρ is a density over an ensemble of many imagined copies, or over your uncertainty about one system. A single trajectory is just one point; Liouville's theorem is a statement about the cloud, not the point.
• "Numerical integrators automatically respect it." Generic Runge–Kutta schemes leak or pump phase volume and drift in energy. You need a symplectic integrator (leapfrog, Verlet) to preserve volume to machine precision over long runs — which is why molecular dynamics insists on them.

Applications

• Statistical mechanics. Justifies the microcanonical and canonical ensembles as stationary distributions, and underwrites the ergodic program connecting time averages to ensemble averages.
• Accelerator and beam physics. Conservation of phase volume is conservation of emittance. Designing a focusing lattice is an exercise in reshaping a fixed-area phase ellipse; cooling techniques (stochastic, electron, laser cooling) are the only ways to genuinely shrink it, because they break the Hamiltonian assumption.
• Plasma and galactic dynamics. The collisionless Vlasov equation is literally Liouville's theorem for a smooth distribution function f, used for tokamak plasmas and for the dynamics of stars in a galaxy.
• Molecular dynamics. Symplectic integrators are chosen precisely so simulations inherit the volume-preserving structure and don't spuriously heat or cool over millions of steps.
• Chaos and mixing theory. Liouville volume conservation plus stretching-and-folding is the engine of chaotic mixing; it sets up the distinction between fine- and coarse-grained entropy and motivates SRB measures.
• Quantum foundations. The von Neumann equation is the structural heir of Liouville's theorem; unitary evolution preserving entropy is the quantum echo of volume preservation, central to debates about thermalization and the black-hole information paradox.

Performance / derivation analysis

The cleanest way to see why the volume is fixed, rather than just verifying it, is to track the infinitesimal volume element under the flow map Φ_t. Volume after time t is

V(t) = ∫ |det J(t)| dV(0),   J = ∂(q(t), p(t)) / ∂(q(0), p(0))

Differentiate the Jacobian determinant. A standard identity gives d/dt ln det J = tr(J⁻¹ J̇) = ∇·v evaluated along the trajectory. We already showed ∇·v = 0 for any Hamiltonian, so:

d/dt (det J) = (∇·v) · det J = 0     ⟹     det J(t) = det J(0) = 1

The Jacobian determinant is pinned at 1 for all time, so V(t) = V(0) — phase volume conserved, full stop. Compare this to the dissipative case H + drag, where ∇·v = −γ and the same identity yields det J(t) = e^(−γt): volume collapses exponentially, which is exactly how trajectories funnel onto an attractor of measure zero.

A self-contained numerical check, integrating an ensemble through a pendulum with a symplectic step:

// Track phase-space area of an ensemble under symplectic (leapfrog) flow.
// Hamiltonian: H = p^2/2 - cos(q)  (pendulum). Hamilton: q' = p, p' = -sin(q).

function leapfrog(q, p, dt) {
  // Symplectic Verlet step -> exactly volume-preserving (det Jacobian = 1)
  const pHalf = p - 0.5 * dt * Math.sin(q);
  const qNew  = q + dt * pHalf;
  const pNew  = pHalf - 0.5 * dt * Math.sin(qNew);
  return [qNew, pNew];
}

// Seed a small parallelogram and follow its two edge vectors.
let q0 = 0.5, p0 = 0.0;
let e1 = [1e-3, 0], e2 = [0, 1e-3];          // edges spanning the blob
const area0 = Math.abs(e1[0] * e2[1] - e1[1] * e2[0]);

const dt = 0.01;
for (let n = 0; n < 200000; n++) {
  // Advance the corner and the two displaced corners, recover edge vectors.
  const c  = leapfrog(q0, p0, dt);
  const c1 = leapfrog(q0 + e1[0], p0 + e1[1], dt);
  const c2 = leapfrog(q0 + e2[0], p0 + e2[1], dt);
  q0 = c[0]; p0 = c[1];
  e1 = [c1[0] - q0, c1[1] - p0];
  e2 = [c2[0] - q0, c2[1] - p0];
}
const area = Math.abs(e1[0] * e2[1] - e1[1] * e2[0]);
console.log(area0.toExponential(6)); // 1.000000e-6
console.log(area.toExponential(6));  // 1.000000e-6  -> conserved to ~1e-12 relative
// The blob is now hugely stretched along the separatrix, yet its AREA is unchanged.

The two edge vectors grow and rotate as the blob smears along the pendulum's separatrix, but their cross product — the area — stays glued to its initial value across 200,000 steps. Swap the symplectic step for a plain forward-Euler step and the same loop reports area drifting by tens of percent: a numerical demonstration that Liouville's theorem is a property of the geometry of the flow, not an accident of any particular integrator.