Part II: Hamiltonian MechanicsTopic 3 of 3

Phase Space & Liouville's Theorem

Phase space is where classical mechanics truly lives. Liouville's theorem guarantees that Hamiltonian flow is incompressible, Poincare recurrence says systems return arbitrarily close to their initial state, and the KAM theorem governs the boundary between order and chaos.

Historical Context

Joseph Liouville proved his theorem in 1838, just a few years after Hamilton's reformulation. The result — that Hamiltonian flow preserves phase-space volume — was recognized by Boltzmann and Gibbs as the foundation of statistical mechanics. Without it, the microcanonical ensemble and the ergodic hypothesis would have no justification.

Henri Poincare proved his recurrence theorem in 1890 while studying the three-body problem. It states that almost every initial condition in a bounded Hamiltonian system will return arbitrarily close to itself. This result led to the famous objection to Boltzmann's H-theorem by Zermelo (1896), and the resolution of this apparent paradox deepened our understanding of statistical mechanics and irreversibility.

The KAM theorem (Kolmogorov 1954, Arnold 1963, Moser 1962) resolved a centuries-old question: does a small perturbation of an integrable system destroy all invariant tori? The answer is no — most tori survive, and chaos is confined to thin layers around resonant tori. This is one of the deepest results in dynamical systems theory.

1. The Structure of Phase Space

For a system with \(n\) degrees of freedom, phase space is the \(2n\)-dimensional space of all possible states \((q_1, \ldots, q_n, p_1, \ldots, p_n)\). Each point represents a complete instantaneous state of the system: both the configuration and the motion.

Hamilton's equations define a flow on phase space: given any initial point \(\mathbf{z}_0 = (q_0, p_0)\) at time \(t_0\), the equations uniquely determine the state at all future (and past) times. The flow map\(\phi_t: \mathbf{z}_0 \mapsto \mathbf{z}(t)\) is the phase-space analog of a velocity field in fluid dynamics.

Key Property: Determinism

By the existence and uniqueness theorem for ODEs, phase-space trajectories never cross. Two distinct initial conditions produce two distinct trajectories for all time. This means phase space has the structure of a foliation by non-intersecting curves.

The Hamiltonian vector field that generates the flow is:

\[X_H = \sum_i \left(\frac{\partial H}{\partial p_i}\frac{\partial}{\partial q_i} - \frac{\partial H}{\partial q_i}\frac{\partial}{\partial p_i}\right)\]

The divergence of this vector field is:

\[\nabla \cdot X_H = \sum_i \left(\frac{\partial^2 H}{\partial q_i \partial p_i} - \frac{\partial^2 H}{\partial p_i \partial q_i}\right) = 0\]

The divergence vanishes identically! This is the basis of Liouville's theorem.

Figure 1. Phase space portrait. Blue ellipses show harmonic oscillator orbits (libration). The gold figure-eight curve is the separatrix of a pendulum, separating bounded libration from unbounded rotation. Arrows indicate the direction of Hamiltonian flow.

2. Liouville's Theorem: The Incompressibility of Hamiltonian Flow

Theorem (Liouville, 1838): The phase-space volume element is preserved under Hamiltonian evolution. If \(\mathcal{V}(t)\) is any region of phase space that evolves under Hamilton's equations, then its volume remains constant:

\[\frac{d}{dt}\int_{\mathcal{V}(t)} d^n q\,d^n p = 0\]

Proof

Consider an infinitesimal phase-space volume element \(d\Gamma = d^n q\,d^n p\). Under Hamiltonian evolution for time \(dt\), the coordinates transform as:

\[q_i \to q_i + \frac{\partial H}{\partial p_i}dt, \quad p_i \to p_i - \frac{\partial H}{\partial q_i}dt\]

The Jacobian of this transformation is:

\[\frac{\partial(q', p')}{\partial(q, p)} = \det(I + J\cdot\nabla^2 H\,dt) = 1 + (\text{tr}\,J\nabla^2 H)\,dt + O(dt^2)\]

The trace vanishes because \(J\) is antisymmetric and \(\nabla^2 H\) is symmetric:\(\text{tr}(J\nabla^2 H) = \sum_i (\partial^2 H/\partial q_i\partial p_i - \partial^2 H/\partial p_i\partial q_i) = 0\). Therefore the Jacobian is \(1\) to all orders in \(dt\), and volumes are preserved. \(\square\)

The Density Function Formulation

Equivalently, if \(\rho(q, p, t)\) is a phase-space density (e.g., an ensemble of systems in statistical mechanics), Liouville's theorem takes the form:

Liouville Equation

\[\frac{\partial \rho}{\partial t} + \{\rho, H\} = 0 \quad \Leftrightarrow \quad \frac{d\rho}{dt} = 0\]

The phase-space density is constant along trajectories — like an incompressible fluid. This is the foundation of the microcanonical ensemble in statistical mechanics.

3. Poincare Recurrence Theorem

Theorem (Poincare, 1890): For a Hamiltonian system confined to a bounded region of phase space, almost every initial condition will return arbitrarily close to itself after a sufficiently long time.

Proof Sketch

Let \(U\) be any open set in the accessible phase space. Consider the sequence of images \(U, \phi_T(U), \phi_{2T}(U), \ldots\) for some fixed \(T\). Since all these sets have the same volume (Liouville) and the total accessible volume is finite, they cannot all be disjoint. So there exist \(n > m \geq 0\) such that\(\phi_{nT}(U) \cap \phi_{mT}(U) \neq \emptyset\).

Applying \(\phi_{-mT}\) to both sides: \(\phi_{(n-m)T}(U) \cap U \neq \emptyset\). So some point in \(U\) returns to \(U\) after time \((n-m)T\). Since \(U\)is arbitrary, almost every point is recurrent. \(\square\)

The Recurrence Paradox

Zermelo (1896) argued that Poincare recurrence contradicts the second law of thermodynamics: if a system always returns to its initial state, how can entropy increase? Boltzmann's resolution: the recurrence time for a macroscopic system is astronomically large — far exceeding the age of the universe. For \(N\)particles, the recurrence time scales as \(e^{cN}\), making it irrelevant for practical thermodynamics.

4. Action-Angle Variables and Canonical Invariants

For an integrable system — one with \(n\) independent conserved quantities in involution — there exists a canonical transformation to action-angle variables \((\theta_i, J_i)\) such that:

\[H = H(J_1, \ldots, J_n)\]

The Hamiltonian depends only on the actions, not the angles. Hamilton's equations then give:

\[\dot{\theta}_i = \frac{\partial H}{\partial J_i} \equiv \omega_i(J) = \text{const}, \quad \dot{J}_i = -\frac{\partial H}{\partial \theta_i} = 0\]

The actions are constants of motion, and the angles increase linearly in time:\(\theta_i(t) = \omega_i t + \theta_i^{(0)}\). Each angle is periodic with period\(2\pi\), so the motion takes place on an \(n\)-dimensional invariant torus \(\mathbb{T}^n\).

The Action Integral

The action variable is defined as:

\[J_i = \frac{1}{2\pi}\oint p_i\,dq_i\]

where the integral is over one complete cycle of the \(i\)-th motion. This is an adiabatic invariant: it remains approximately constant even when external parameters are changed slowly. This is the classical precursor of the quantum number (Bohr-Sommerfeld quantization: \(J_i = n_i \hbar\)).

Example: Harmonic Oscillator

\(H = p^2/(2m) + m\omega^2 q^2/2\). The orbit in phase space is an ellipse with area \(2\pi J = 2\pi E/\omega\). So \(J = E/\omega\) and \(H = \omega J\).

The angle variable is \(\theta = \omega t\), uniformly increasing. The Bohr-Sommerfeld condition \(J = n\hbar\) gives \(E_n = n\hbar\omega\) — the correct quantum energy levels (up to the zero-point energy).

5. The KAM Theorem and the Onset of Chaos

What happens when an integrable system is perturbed? If \(H = H_0(J) + \varepsilon H_1(\theta, J)\), do the invariant tori of \(H_0\) survive?

KAM Theorem (Kolmogorov-Arnold-Moser): If the perturbation \(\varepsilon\)is sufficiently small and the unperturbed system is non-degenerate (\(\det(\partial^2 H_0/\partial J_i \partial J_j) \neq 0\)), then most invariant tori survive, being only slightly deformed. The tori that survive are those whose frequency ratios \(\omega_i/\omega_j\) are sufficiently irrational (satisfying a Diophantine condition).

The tori that are destroyed are the resonant ones, where\(\sum_i n_i \omega_i = 0\) for integer \(n_i\). Near these destroyed tori, chaotic motion appears. For small \(\varepsilon\), the chaotic regions are thin layers between surviving KAM tori. As \(\varepsilon\) increases, more tori break, the chaotic layers widen, and eventually overlap — leading to large-scale (global) chaos.

The Standard Map

The paradigmatic model for studying the KAM transition is the standard map (Chirikov, 1969):

\[p_{n+1} = p_n + K\sin(\theta_n), \quad \theta_{n+1} = \theta_n + p_{n+1}\]

For \(K = 0\), all orbits are on invariant tori (straight lines in the\((\theta, p)\) plane). As \(K\) increases, islands of regular motion and seas of chaos appear. At \(K \approx 0.97\) (the "critical" value), the last spanning KAM torus (the "golden ratio" torus) breaks, and global diffusion becomes possible.

6. Chaos in Hamiltonian Systems

Hamiltonian chaos is subtly different from dissipative chaos (like the Lorenz attractor). By Liouville's theorem, Hamiltonian systems cannot have attractors — phase-space volumes are preserved. Instead, chaotic motion fills regions of phase space ergodically.

The hallmark of chaos is sensitive dependence on initial conditions, quantified by Lyapunov exponents. For a Hamiltonian system with \(n\) degrees of freedom, there are \(2n\) Lyapunov exponents\(\lambda_1 \geq \lambda_2 \geq \ldots \geq \lambda_{2n}\). By Liouville's theorem, their sum vanishes:

\[\sum_{i=1}^{2n} \lambda_i = 0\]

Furthermore, they come in pairs: \(\lambda_i = -\lambda_{2n+1-i}\). If the largest Lyapunov exponent \(\lambda_1 > 0\), nearby trajectories diverge exponentially:\(|\delta\mathbf{z}(t)| \sim |\delta\mathbf{z}(0)|\,e^{\lambda_1 t}\).

For integrable systems, all Lyapunov exponents are zero (power-law divergence at most). The presence of positive Lyapunov exponents is the definitive signature of chaos.

Python Simulation: Liouville, KAM, and Chaos

We demonstrate Liouville's theorem by tracking a cloud of phase-space points, visualize the KAM transition using the standard map, compute Poincare sections for the Henon-Heiles system, and measure Lyapunov exponents.

Phase Space: Liouville Theorem, KAM Transition, and Chaos

Python

script.py232 lines

import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp

# ============================================================
# Part 1: Liouville's theorem - phase space volume conservation
# ============================================================
fig, axes = plt.subplots(1, 3, figsize=(18, 5))

# Harmonic oscillator: H = (p^2 + q^2)/2
def ho(t, z):
    q, p = z[::2], z[1::2]
    dq = p
    dp = -q
    out = np.empty_like(z)
    out[::2] = dq
    out[1::2] = dp
    return out

# Initial cloud: circle in phase space
N_pts = 500
theta_init = np.linspace(0, 2*np.pi, N_pts, endpoint=False)
r_cloud = 0.3
q0_cloud = 1.0 + r_cloud * np.cos(theta_init)
p0_cloud = 0.5 + r_cloud * np.sin(theta_init)

z0 = np.empty(2 * N_pts)
z0[::2] = q0_cloud
z0[1::2] = p0_cloud

times_plot = [0, 1.0, 2.5]
for idx, t_end in enumerate(times_plot):
    if t_end == 0:
        q_pts = q0_cloud
        p_pts = p0_cloud
    else:
        sol = solve_ivp(ho, [0, t_end], z0, max_step=0.01, rtol=1e-10)
        q_pts = sol.y[::2, -1]
        p_pts = sol.y[1::2, -1]

ax = axes[idx]
    ax.fill(q_pts, p_pts, alpha=0.3, color='royalblue')
    ax.plot(q_pts, p_pts, 'b.', markersize=1)
    ax.set_xlabel('q')
    ax.set_ylabel('p')
    ax.set_title('t = {:.1f}'.format(t_end))
    ax.set_xlim(-2.5, 2.5)
    ax.set_ylim(-2.5, 2.5)
    ax.set_aspect('equal')
    ax.grid(True, alpha=0.3)

# Compute area using shoelace formula
    area = 0.5 * np.abs(np.sum(q_pts * np.roll(p_pts, -1) - np.roll(q_pts, -1) * p_pts))
    ax.set_title('t = {:.1f}, Area = {:.4f}'.format(t_end, area))

plt.suptitle("Liouville's Theorem: Phase Space Area Conservation", fontsize=14, y=1.02)
plt.tight_layout()
plt.savefig('output.png', dpi=150, bbox_inches='tight')
plt.close()

print("--- Liouville's Theorem ---")
print("Initial area of phase-space cloud: {:.6f}".format(np.pi * r_cloud**2))
print("Area is preserved at all times (shape changes, volume doesn't)")

# ============================================================
# Part 2: Standard Map - KAM transition
# ============================================================
fig2, axes2 = plt.subplots(1, 3, figsize=(18, 5))

K_values = [0.2, 0.9, 2.0]

for idx, K in enumerate(K_values):
    ax = axes2[idx]

for p0 in np.linspace(0, 2*np.pi, 30):
        for th0 in [0.1, 1.0, 2.0, 3.0, 4.0, 5.0]:
            theta, p = th0, p0
            thetas, ps = [], []
            for _ in range(500):
                p = p + K * np.sin(theta)
                theta = theta + p
                # Fold into [0, 2pi)
                theta = theta % (2*np.pi)
                p = p % (2*np.pi)
                thetas.append(theta)
                ps.append(p)
            ax.plot(thetas, ps, 'b.', markersize=0.15, alpha=0.5)

ax.set_xlabel('theta')
    ax.set_ylabel('p')
    ax.set_title('Standard Map, K = {:.1f}'.format(K))
    ax.set_xlim(0, 2*np.pi)
    ax.set_ylim(0, 2*np.pi)
    ax.grid(True, alpha=0.2)

plt.suptitle('KAM Theorem: Destruction of Invariant Tori', fontsize=14, y=1.02)
plt.tight_layout()
plt.savefig('output2.png', dpi=150, bbox_inches='tight')
plt.close()

print("\n--- Standard Map (KAM Transition) ---")
print("K = 0.2: Most invariant tori survive (regular motion)")
print("K = 0.9: Islands of regularity in seas of chaos")
print("K = 2.0: Large-scale chaos, few surviving tori")

# ============================================================
# Part 3: Henon-Heiles - order to chaos transition
# ============================================================
def henon_heiles(t, z):
    x, y, px, py = z
    dxdt = px
    dydt = py
    dpx = -x - 2*x*y
    dpy = -y - x**2 + y**2
    return [dxdt, dydt, dpx, dpy]

fig3, axes3 = plt.subplots(1, 3, figsize=(18, 5))
energies_hh = [1.0/20, 1.0/10, 1.0/6.5]
titles_hh = ['E=0.05 (regular)', 'E=0.10 (mixed)', 'E=0.154 (chaotic)']

for idx, (E_hh, title_hh) in enumerate(zip(energies_hh, titles_hh)):
    ax = axes3[idx]
    x_sec_list, px_sec_list = [], []

for x0 in np.linspace(-0.5, 0.5, 40):
        for px0 in np.linspace(-0.5, 0.5, 40):
            py2 = 2*E_hh - px0**2 - x0**2
            if py2 < 0:
                continue
            py0 = np.sqrt(py2)

try:
                sol = solve_ivp(henon_heiles, [0, 300], [x0, 0, px0, py0],
                               max_step=0.1, rtol=1e-9, atol=1e-11)
            except:
                continue

y_arr = sol.y[1]
            py_arr = sol.y[3]
            for j in range(1, len(y_arr)):
                if y_arr[j-1] < 0 and y_arr[j] >= 0 and py_arr[j] > 0:
                    frac = -y_arr[j-1] / (y_arr[j] - y_arr[j-1])
                    x_cross = sol.y[0][j-1] + frac * (sol.y[0][j] - sol.y[0][j-1])
                    px_cross = sol.y[2][j-1] + frac * (sol.y[2][j] - sol.y[2][j-1])
                    x_sec_list.append(x_cross)
                    px_sec_list.append(px_cross)

ax.plot(x_sec_list, px_sec_list, 'b.', markersize=0.2, alpha=0.5)
    ax.set_xlabel('x')
    ax.set_ylabel('px')
    ax.set_title(title_hh)
    ax.grid(True, alpha=0.2)

plt.suptitle('Henon-Heiles: Transition to Chaos', fontsize=14, y=1.02)
plt.tight_layout()
plt.savefig('output3.png', dpi=150, bbox_inches='tight')
plt.close()

# ============================================================
# Part 4: Lyapunov exponent computation
# ============================================================
print("\n--- Lyapunov Exponents ---")

def compute_lyapunov_ho(T_final=100):
    """Harmonic oscillator (integrable) -> lambda = 0"""
    dt = 0.01
    N = int(T_final / dt)
    z = np.array([1.0, 0.5])
    dz = np.array([1e-8, 0.0])

lyap_sum = 0
    for i in range(N):
        # Evolve both
        q, p = z
        z_new = np.array([q*np.cos(dt) + p*np.sin(dt),
                          -q*np.sin(dt) + p*np.cos(dt)])
        dq, dp = z + dz
        dz_new = np.array([dq*np.cos(dt) + dp*np.sin(dt),
                           -dq*np.sin(dt) + dp*np.cos(dt)]) - z_new

norm_dz = np.linalg.norm(dz_new)
        lyap_sum += np.log(norm_dz / np.linalg.norm(dz))
        dz = dz_new / norm_dz * np.linalg.norm(dz)
        z = z_new

return lyap_sum / (N * dt)

def compute_lyapunov_sm(K, T_final=10000):
    """Standard map Lyapunov exponent."""
    theta, p = 0.1, 0.5
    dtheta, dp_pert = 1e-8, 0.0

lyap_sum = 0
    for _ in range(T_final):
        # Evolve reference
        p_new = p + K * np.sin(theta)
        theta_new = theta + p_new

# Evolve perturbation (tangent map)
        dp_new = dp_pert + K * np.cos(theta) * dtheta
        dtheta_new = dtheta + dp_new

norm = np.sqrt(dtheta_new**2 + dp_new**2)
        lyap_sum += np.log(norm / np.sqrt(dtheta**2 + dp_pert**2 + 1e-30))

# Renormalize
        scale = np.sqrt(dtheta**2 + dp_pert**2)
        dtheta = dtheta_new / norm * scale
        dp_pert = dp_new / norm * scale
        theta = theta_new % (2*np.pi)
        p = p_new % (2*np.pi)

return lyap_sum / T_final

lam_ho = compute_lyapunov_ho()
print("Harmonic oscillator (integrable): lambda = {:.6f} (should be ~0)".format(lam_ho))

for K in [0.2, 0.9, 2.0, 5.0]:
    lam_sm = compute_lyapunov_sm(K)
    print("Standard map K={:.1f}: lambda = {:.4f}".format(K, lam_sm))

print("\nPositive Lyapunov exponent = chaos (exponential divergence)")
print("Zero Lyapunov exponent = regular motion (power-law divergence)")

print("\nSimulation complete. Three plots generated.")
print("Plot 1: Liouville theorem - phase space area conservation")
print("Plot 2: Standard map - KAM transition at K=0.2, 0.9, 2.0")
print("Plot 3: Henon-Heiles Poincare sections at three energies")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

Summary

Phase space is the \(2n\)-dimensional space of \((q, p)\) — the natural arena for Hamiltonian mechanics.
Liouville's theorem: Hamiltonian flow preserves phase-space volume (incompressible flow).
Poincare recurrence: bounded Hamiltonian systems return arbitrarily close to initial conditions.
Action-angle variables linearize integrable systems; actions are adiabatic invariants and precursors to quantum numbers.
KAM theorem: most invariant tori survive small perturbations; chaos appears near resonances.
Lyapunov exponents quantify chaos; they come in pairs summing to zero (Liouville constraint).
The standard map and Henon-Heiles system are paradigmatic models for studying the transition from order to chaos.

← Hamilton's Equations Back to Course Overview →