Basic statistical mechanics

After deriving the equations of motion for a system of \(N\) classical particles, we would next like to solve for their trajectories and use them to evaluate thermodynamic properties of a macroscopic system, such as the temperature, energy, pressure, etc. If we can connect an observable macroscopic property \(A\) to some function of the microstate of the system, \(a(\textbf{x})\), we might associate the measured value of \(A\) with a time average, denoted by \(\bar{a}_{\tau}\) and defined as

(22)\[ A = \bar{a}_{\tau} = \frac{1}{\tau} \int_0^{\tau} a(\textbf{x}(t)) \mathrm{d} t, \]

where \(\tau\) is some long time scale compared to the molecular or atomic motion.

But now we find ourselves at an impasse. It is entirely clear what rules determine the evolution of the system in time, but solving the equations analytically is impossible. A numerical solution is similarly unfeasible for \(\sim 10^{23}\) particles, as we discussed in the previous chapter. The core idea of statistical mechanics is to provide an alternative approach to obtain macroscopic thermodynamic observables.

Ensembles

It was J.W. Gibbs that realized that thermodynamic macroscopic properties are not very sensitive to the details of the specific microscopic state of the system. Consider for example the temperature. You have learned in a preceding statistical thermodynamics course that each momentum degree of freedom contributes \(\frac{1}{2} k_B T \) to the average thermal kinetic energy (the equipartition theorem). Thus, the temperature can be obtained from the average kinetic energy of the system. But there are many ways to choose the momenta of the particles such that the average kinetic energy is the same. In other words, the temperature is insensitive to the specific microstate of the system, as long as it gives the same average kinetic energy.

Following this reasoning, and to make progress in connecting microscopic properties with thermodynamic observables, Gibbs introduced the concept of an ensemble. An ensemble is a collection of \(\mathcal{N}\) systems that obey the same microscopic interactions and are characterized by a common set of three thermodynamic properties, such as the number of particles \(N\), the volume \(V\), the total energy \(E\), the temperature \(T\), the pressure \(P\) or the chemical potential \(\mu\). Each system in the ensemble evolves in time from different initial conditions and is characterized by a unique microscopic state defined by its path in phase space, \(\textbf{x}(t)\), defined in Eq. (7). The entire ensemble can then be represented as a “cloud” of points in phase space.

Fig. 3 An ensemble: \(\mathcal{N}\) copies of the same system, with the same \(N\), \(V\) and \(E\), each in its own microstate \(\textbf{x}^{(k)}(t)\).

We define the function \(f(\textbf{x},t)\), called the phase space distribution function, as the probability density of observing a system in the ensemble in a specific microstate at time \(t\). Since it is a probability density, it has to fulfill the properties:

(23)\[\begin{align} \int f(\textbf{x},t) \, \mathrm{d} \textbf{x} &= 1, & f(\textbf{x},t) & \ge 0. \end{align}\]

In other words, \(f(\textbf{x},t) \, \mathrm{d} \textbf{x}\) represents the probability that a system from the ensemble is in the volume element between \(\textbf{x}\) and \(\textbf{x}+\mathrm{d}\textbf{x}\) at time \(t\) (or, alternatively, the fraction of systems therein).

Gibbs then postulated that instead of using the time average of Eq. (22), macroscopic properties can be obtained by freezing the ensemble at some time instance and taking the average of the desired observable over all the systems, called an ensemble average. More formally, the ensemble average, denoted as \(\langle a \rangle_f\), is defined as

(24)\[ A(t) = \langle a \rangle_f = \int f(\textbf{x},t) a(\textbf{x}) \, \mathrm{d} \textbf{x}. \]

Note the difference between Eqs. (22) and (24), the former is an integral over time while the latter is an integral over phase space. This distinction is important enough that we have used a bar to denote the time average but brackets to denote the ensemble average. Before moving forward, and discussing what we can do with ensemble averages, we will briefly investigate an important property of the phase space distribution function, described by the Liouville equation.

Active learning

For a single particle of mass \(m\) in a harmonic trap of frequency \(\omega_0\), calculate the equilibrium average of the total energy in the canonical ensemble.

Hint: In that case,

\[f(q,p) = \frac{e^{-\beta \mathcal{H}(q,p)}}{\mathcal{Z}} \]

You will need the Gaussian integrals, for \(a > 0\),

\[ \int_{-\infty}^{\infty} e^{-a x^2} \, \mathrm{d}x = \sqrt{\frac{\pi}{a}}, \qquad \int_{-\infty}^{\infty} x^2 e^{-a x^2} \, \mathrm{d}x = \frac{1}{2a} \sqrt{\frac{\pi}{a}}. \]

Solution

The ensemble average is

\[ \langle \mathcal{H} \rangle = \frac{1}{\mathcal{Z}} \int \mathcal{H}(q,p) \, e^{-\beta \mathcal{H}(q,p)} \, \mathrm{d}q \, \mathrm{d}p, \qquad \mathcal{Z} = \int e^{-\beta \mathcal{H}(q,p)} \, \mathrm{d}q \, \mathrm{d}p = \frac{2\pi}{\beta \omega_0} \]

The Hamiltonian is

\[ \mathcal{H} = \frac{p^2}{2m} + \frac{1}{2} m \omega_0^2 q^2 \]

After doing some Gaussian integrals, we get:

\[ \langle \mathcal{H} \rangle = \beta^{-1} \]

Liouville equation

For clarity, we will prove the Liouville equation for a single particle, for which \( f(\textbf{x},t) = f(q,p,t) \), but it can be generalized to \(N\) particles. The theorem states that the phase space distribution function is conserved along the trajectory, i.e., that

(25)\[ \frac{\mathrm{d} f}{\mathrm{d} t} = \frac{\partial f}{\partial t} + \frac{\partial f }{\partial q} \dot{q} + \frac{\partial f }{\partial p} \dot{p} = 0. \]

The proof has two steps. The first uses the equations of motion, the second only counts systems.

Step 1: a property of Hamilton’s equations

Hamilton’s equations assign to every point \(\textbf{x}=[q,p]\) of phase space a velocity \(\dot{\textbf{x}}=[\dot{q},\dot{p}]\), which tells us where a system at that point goes next. This velocity has a special property.

Active learning

Show that

\[ \frac{\partial \dot{q}}{\partial q} + \frac{\partial \dot{p}}{\partial p} = 0. \]

Hint: Use Hamilton’s equations.

Solution

Using Hamilton’s equations, Eq. (13), \(\dot{q} = \partial \mathcal{H} / \partial p\) and \(\dot{p} = - \partial \mathcal{H} / \partial q\), so

(26)\[ \frac{\partial \dot{q}}{\partial q} + \frac{\partial \dot{p}}{\partial p} = \frac{ \partial ^2 \mathcal{H}}{\partial q \partial p} - \frac{ \partial ^2 \mathcal{H}}{\partial p \partial q} = 0. \]

In other words, the divergence of the velocity in phase space vanishes, \(\nabla \cdot \dot{\textbf{x}} = 0\). This is the only place in the proof where the specific form of the equations of motion enters. Everything that follows holds for any dynamics with this property.

Step 2: counting systems

We now consider some volume in phase space \(\Omega\) and define the number of systems in the ensemble whose microstate is in this volume as

\[ \mathcal{N}_{\Omega} = \mathcal{N} \int_{\Omega} f(\textbf{x},t) \, \mathrm{d} \textbf{x} = \mathcal{N} \int_{\Omega} f(q,p,t) \, \mathrm{d} q \, \mathrm{d} p. \]

The rate of change of the number of systems in \(\Omega\) is given by

(27)\[ \frac{\mathrm{d} \mathcal{N}_{\Omega}}{\mathrm{d} t} = \mathcal{N} \int_{\Omega} \frac{\partial f}{\partial t} \, \mathrm{d} q \, \mathrm{d} p. \]

But, equivalently, since systems are neither created nor destroyed, the number in \(\Omega\) changes only when systems cross its boundary.

The flux of systems leaving the volume through a surface element \(\mathrm{d} \mathcal{A}\) can be integrated to give the rate,

(28)\[ \frac{\mathrm{d} \mathcal{N}_{\Omega}}{\mathrm{d} t} = - \mathcal{N} \oint_{\mathcal{A}} f \dot{\textbf{x}} \cdot \textbf{n} \, \mathrm{d} \mathcal{A}, \]

where \(\textbf{n}\) is the unit vector normal to the surface element, pointing outward. For a single particle, phase space is the \((q,p)\) plane, \(\Omega\) is a region in it, \(\mathcal{A}\) is the closed curve bounding it and \(\mathrm{d} \mathcal{A}\) is a line element of that curve. Only the component of the velocity \(\dot{\textbf{x}} = [\dot{q},\dot{p}]\) along \(\textbf{n}\) carries systems out of \(\Omega\). For example, for the small rectangle in Fig. 4, the right side contributes \(f \dot{q} \, \mathrm{d}p\) to the integral and the top side \(f \dot{p} \, \mathrm{d}q\).

Fig. 4 A small rectangle \(\Omega\) in phase space with sides \(\mathrm{d}q\) and \(\mathrm{d}p\). The flux of systems through its right side, where the outward normal \(\textbf{n}\) is along \(q\), is \(f \dot{q} \, \mathrm{d}p\), and through its top side \(f \dot{p} \, \mathrm{d}q\).

We can now use the divergence theorem, which turns the integral over the boundary curve into an integral over the region,

(29)\[ \frac{\mathrm{d} \mathcal{N}_{\Omega}}{\mathrm{d} t} = - \mathcal{N} \int_{\Omega} \nabla \cdot ( f \dot{\textbf{x}} ) \, \mathrm{d} \Omega = - \mathcal{N} \int_{\Omega} \left[ \frac{\partial (f \dot{q})}{\partial q} + \frac{\partial (f \dot{p})}{\partial p} \right] \, \mathrm{d} q \, \mathrm{d} p. \]

Eq.(29) and Eq. (27) should be equal for every volume element \(\mathrm{d} q \, \mathrm{d} p \), giving

(30)\[ \frac{\partial f}{\partial t} + \left[ \frac{\partial (f \dot{q})}{\partial q} + \frac{\partial (f \dot{p})}{\partial p} \right] = 0. \]

Active learning

Show that the term in square brackets can be written as

\[ \frac{\partial (f \dot{q})}{\partial q} + \frac{\partial (f \dot{p})}{\partial p} = \frac{\partial f }{\partial q} \dot{q} + \frac{\partial f }{\partial p} \dot{p}. \]

Solution

The term in square brackets can be rearranged

\[ \frac{\partial (f \dot{q})}{\partial q} + \frac{\partial (f \dot{p})}{\partial p} = f \left( \frac{\partial \dot{q}}{\partial q} + \frac{\partial \dot{p}}{\partial p} \right) + \frac{\partial f }{\partial q} \dot{q} + \frac{\partial f }{\partial p} \dot{p}, \]

and the first term is zero by Eq. (26).

Based on the solution of the above exercise, Eq. (30) reduces to Eq. (25), which completes the proof of the Liouville equation.

Using the definition of the Poisson bracket between \(f\) and the Hamiltonian, Eq. (16), we can write

\[ \frac{\partial f }{\partial q} \dot{q} + \frac{\partial f }{\partial p} \dot{p} = \frac{\partial f }{\partial q} \frac{\partial \mathcal{H}}{\partial p} - \frac{\partial f }{\partial p} \frac{\partial \mathcal{H}}{\partial q} = \{ f, \mathcal{H} \}, \]

which leads to another useful form of the Liouville equation,

(31)\[ \frac{\partial f}{\partial t} + \{ f, \mathcal{H} \} = 0. \]

What the Liouville equation says

Take one system of the ensemble and follow it along its trajectory \(\textbf{x}(t)\). At every instant, evaluate the phase space distribution at the point where this system is, \(f(\textbf{x}(t),t)\). The derivative in Eq. (25) is the rate of change of this quantity: the first term accounts for the explicit dependence of \(f\) on time, the other two for the motion of the system. The Liouville equation says that it is conserved along the dynamics of each system,

\[ f(\textbf{x}(t),t) = f(\textbf{x}(0),0). \]

This does not mean that \(f\) does not change in time. The distribution as a whole moves in phase space, as Fig. 5 shows, and in general also changes its shape, as Fig. 6 shows.

Fig. 5 An ensemble of harmonic oscillators (\(m = \omega_0 = 1\)) whose distribution \(f\) (color) starts in a square in phase space, peaked towards one corner. A few systems of the ensemble are shown as points, four of them with their trajectories. The distribution rotates, but \(f\) along each trajectory does not change.

Fig. 6 The same for an ensemble of pendulums, \(\mathcal{H} = \frac{p_\phi^2}{2 m l^2} + m g l (1 - \cos{\phi})\) as in the previous chapter, with \(m = l = g = 1\). The period depends on the energy, so the distribution changes its shape, but not its area. The purple system starts with the highest energy and ends up in the tail.

Solving the Liouville equation

Can we solve Eq. (25) for \(f(\textbf{x},t)\)? It is a first-order linear partial differential equation, and the curves along which its solution is constant, its characteristics, are precisely the trajectories of the systems. Solving it is therefore equivalent to solving the equations of motion: if we know where every system came from, \(f(\textbf{x},t) = f(\textbf{x}(0),0)\), which is how the animations above were drawn. But then we gain nothing beyond the trajectories themselves. Instead, the Liouville equation is useful in two other ways (and more that we will not discuss here).

1. Equilibrium. We can find distributions that do not change in time at all, \(\partial f / \partial t = 0\), without solving the equations of motion. This is the main use of the Liouville equation in this course, and the subject of the next section.

2. Moments. Multiplying Eq. (25) by \(q\), \(p\), \(q^2\) and so on, and integrating over phase space, we obtain ordinary differential equations for the averages \(\langle q \rangle\), \(\langle p \rangle\), \(\langle q^2 \rangle\) and so on, without ever obtaining the full distribution. For the harmonic oscillator, for example, \(\mathrm{d} \langle q \rangle / \mathrm{d}t = \langle p \rangle / m\) and \(\mathrm{d} \langle p \rangle / \mathrm{d}t = - m \omega_0^2 \langle q \rangle\), so the average position and momentum of the ensemble obey Newton’s equations exactly.

Phase space incompressibility

A second consequence of the property of Hamilton’s equations shown in Step 1 above, Eq. (26), concerns volumes in phase space rather than densities. We denote an initial point in phase space as \(\textbf{x}_0 = [q_0,p_0]\). The point in phase space that results from the dynamics starting from \(\textbf{x}_0\) after time \(t\) is \(\textbf{x}_t = [q_t,p_t]\). In the previous chapter, we learned that the trajectory is uniquely defined by the initial conditions, meaning that \(\textbf{x}_t\) can be written as a function of the initial conditions,

\[ \textbf{x}_t = g(\textbf{x}_0). \]

Thus, we can consider the evolution of the trajectory as a coordinate transformation. Assuming we know the function \(g\), we need to evaluate how the volume elements change when we transform from \(\mathrm{d} \textbf{x}_0\) to \(\mathrm{d} \textbf{x}_t\). J. Liouville proved that

(32)\[ \mathrm{d} \textbf{x}_0 = \mathrm{d} \textbf{x}_t, \]

i.e., that the phase space volume element is conserved along the trajectory. Its shape might change, but not the volume. In Fig. 5, the square of initial conditions rotates, but its area is conserved.

Why is it important? We will use it in the next section. There, we will postulate that at equilibrium all microstates of an isolated system are equally probable, i.e., that \(f\) is uniform on the surface of constant energy. For this to be consistent with the dynamics, the equations of motion must not squeeze or stretch phase space: a distribution that is uniform at \(t=0\) must remain uniform at later times. Eq. (32) guarantees exactly that. It also tells us that the phase space volume \(\mathrm{d} \textbf{x}\) is the natural measure for counting microstates.

To prove Eq. (32), again for a single particle, we remember that two volume elements of a coordinate transformation are connected through

(33)\[ \mathrm{d} \textbf{x}_t = \mathcal{J}(\textbf{x}_t; \textbf{x}_0) \mathrm{d} \textbf{x}_0, \]

where the Jacobian \(\mathcal{J}\) is the determinant of the matrix of partial derivatives,

\[\begin{split} \mathcal{J}(\textbf{x}_t; \textbf{x}_0) = \det{J}, \qquad J= \begin{pmatrix} \frac{\partial q_t}{\partial q_0} & \frac{\partial q_t}{\partial p_0} \\ \frac{\partial p_t}{\partial q_0} & \frac{\partial p_t}{\partial p_0} \end{pmatrix}. \end{split}\]

At \(t=0\) the transformation is the identity, so \(J\) is the unit matrix and \(\mathcal{J}=1\). We will show that \(\mathcal{J}\) does not change in time. The time derivative of the first element of \(J\) is

\[ \frac{\mathrm{d}}{\mathrm{d}t} \frac{\partial q_t}{\partial q_0} = \frac{\partial \dot{q}_t}{\partial q_0} = \frac{\partial \dot{q}_t}{\partial q_t} \frac{\partial q_t}{\partial q_0} + \frac{\partial \dot{q}_t}{\partial p_t} \frac{\partial p_t}{\partial q_0}, \]

where we used the chain rule, since \(\dot{q}_t = \partial \mathcal{H} / \partial p\) is a function of \(\textbf{x}_t\). Doing the same for the other three elements, we find that the whole matrix obeys

\[\begin{split} \frac{\mathrm{d} J}{\mathrm{d}t} = A J, \qquad A= \begin{pmatrix} \frac{\partial \dot{q}_t}{\partial q_t} & \frac{\partial \dot{q}_t}{\partial p_t} \\ \frac{\partial \dot{p}_t}{\partial q_t} & \frac{\partial \dot{p}_t}{\partial p_t} \end{pmatrix}. \end{split}\]

Active learning

Show that for a \(2 \times 2\) matrix obeying \(\frac{\mathrm{d} J}{\mathrm{d}t} = A J\), the determinant obeys

\[ \frac{\mathrm{d}}{\mathrm{d}t} \det{J} = \left( A_{11} + A_{22} \right) \det{J}. \]

Solution

Writing \(\det{J} = J_{11} J_{22} - J_{12} J_{21}\) and differentiating,

\[ \frac{\mathrm{d}}{\mathrm{d}t} \det{J} = \dot{J}_{11} J_{22} + J_{11} \dot{J}_{22} - \dot{J}_{12} J_{21} - J_{12} \dot{J}_{21}. \]

The four entries of \(\dot{J} = A J\) are

(34)\[\begin{align} \dot{J}_{11} &= A_{11} J_{11} + A_{12} J_{21}, & \dot{J}_{12} &= A_{11} J_{12} + A_{12} J_{22}, \\ \dot{J}_{21} &= A_{21} J_{11} + A_{22} J_{21}, & \dot{J}_{22} &= A_{21} J_{12} + A_{22} J_{22}. \end{align}\]

Substituting them,

\[ \frac{\mathrm{d}}{\mathrm{d}t} \det{J} = \left( A_{11} J_{11} + A_{12} J_{21} \right) J_{22} + J_{11} \left( A_{21} J_{12} + A_{22} J_{22} \right) - \left( A_{11} J_{12} + A_{12} J_{22} \right) J_{21} - J_{12} \left( A_{21} J_{11} + A_{22} J_{21} \right). \]

Collecting the terms by the elements of \(A\),

\[ \frac{\mathrm{d}}{\mathrm{d}t} \det{J} = A_{11} \left( J_{11} J_{22} - J_{12} J_{21} \right) + A_{22} \left( J_{11} J_{22} - J_{12} J_{21} \right) + A_{12} \left( J_{21} J_{22} - J_{22} J_{21} \right) + A_{21} \left( J_{11} J_{12} - J_{12} J_{11} \right). \]

The last two brackets vanish, and the first two are both \(\det{J}\), so

\[ \frac{\mathrm{d}}{\mathrm{d}t} \det{J} = \left( A_{11} + A_{22} \right) \det{J}. \]

The sum \(A_{11} + A_{22}\) is exactly the quantity of Step 1, Eq. (26), evaluated at \(\textbf{x}_t\). Therefore,

\[ \frac{\mathrm{d}}{\mathrm{d} t} \mathcal{J}(\textbf{x}_t;\textbf{x}_0) = \left( \frac{\partial \dot{q}_t}{\partial q_t} + \frac{\partial \dot{p}_t}{\partial p_t} \right) \mathcal{J}(\textbf{x}_t;\textbf{x}_0) = 0, \]

and since \(\mathcal{J}=1\) initially, it stays equal to one along the trajectory, thus proving Eq. (32). For \(N\) particles, \(J\) is a \(6N \times 6N\) matrix and the same result follows from the general identity \(\frac{\mathrm{d}}{\mathrm{d}t} \det{J} = \mathrm{Tr}[A] \det{J}\), where the trace is the sum of \(\frac{\partial \dot{q}_i}{\partial q_i} + \frac{\partial \dot{p}_i}{\partial p_i}\) over all degrees of freedom, each pair cancelling as in Step 1.

The Liouville equation and phase space incompressibility are two sides of the same coin. Consider a small volume element that moves along with the systems in it. Since systems are neither created nor destroyed, and since trajectories never cross, the number of systems in it, \(f \, \mathrm{d} \textbf{x}\), is constant. The Liouville equation says that the first factor is constant along the trajectory, incompressibility says that the second one is. Given the first, the second follows, and vice versa.

If you didn’t follow the derivation, do not worry! You will now demonstrate phase space incompressibility for a simple example.

Active learning

Show that the Eq. (32) is correct for the case of a single particle of mass \(m\) in a harmonic trap of frequency \(\omega_0\). In this case, \(\textbf{x} = [q,p]\). You can use the results of the previous chapter,

\[ q(t) = q_0 \cos{\omega_0 t} + \frac{p_0}{m \omega_0} \sin{\omega_0 t} \]

\[ p(t) = p_0 \cos{\omega_0 t} - m \omega_0 q_0 \sin{\omega_0 t} \]

Hint: Compute the Jacobian of the coordinate transformation and show that it is unity.

Solution

In the case of a single particle, the Jacobian is a 2x2 matrix,

\[\begin{split} J= \begin{pmatrix} \frac{\partial q_t}{\partial q_0} & \frac{\partial q_t}{\partial p_0} \\ \frac{\partial p_t}{\partial q_0} & \frac{\partial p_t}{\partial p_0} \end{pmatrix} = \begin{pmatrix} \cos{\omega_0 t} & \sin(\omega_0 t) / m \omega_0 \\ - m \omega_0 \sin(\omega_0 t ) & \cos{\omega_0 t} \end{pmatrix} \end{split}\]

whose determinant is \(\det{J} = \cos^2{\omega_0 t} + \sin^2(\omega_0 t ) = 1\).

Extra reading: a damped oscillator

Does volume conservation hold for any Hamiltonian, even one that depends on time? Consider a damped harmonic oscillator, with the equation of motion

\[ \ddot{q} = - \gamma \dot{q} - \omega_0^2 q, \]

where \(\gamma\) is a friction coefficient. Friction is not a conservative force, but the equation of motion can still be obtained from a Lagrangian,

\[ \mathcal{L}(q,\dot{q},t) = e^{\gamma t} \left( \frac{1}{2} m \dot{q}^2 - \frac{1}{2} m \omega_0^2 q^2 \right), \]

which depends explicitly on time. Indeed, the Euler-Lagrange equation, Eq. (10), gives

\[ \frac{\mathrm{d}}{\mathrm{d}t} \left( m e^{\gamma t} \dot{q} \right) = m e^{\gamma t} \left( \ddot{q} + \gamma \dot{q} \right) = - m \omega_0^2 e^{\gamma t} q, \]

and dividing by \(m e^{\gamma t}\) recovers the equation of motion. The momentum conjugate to \(q\) is

\[ p = \frac{\partial \mathcal{L}}{\partial \dot{q}} = m e^{\gamma t} \dot{q}, \]

which is not the linear momentum \(m \dot{q}\). The Hamiltonian, Eq. (12), is

\[ \mathcal{H}(q,p,t) = p \dot{q} - \mathcal{L} = e^{-\gamma t} \frac{p^2}{2m} + e^{\gamma t} \frac{1}{2} m \omega_0^2 q^2. \]

Unlike the cases of the previous chapter, it depends explicitly on time and it is not \(K + U\). Hamilton’s equations, Eq. (13), are

\[ \dot{q} = \frac{\partial \mathcal{H}}{\partial p} = e^{-\gamma t} \frac{p}{m}, \qquad \dot{p} = - \frac{\partial \mathcal{H}}{\partial q} = - e^{\gamma t} m \omega_0^2 q, \]

and differentiating the first and substituting the second recovers the equation of motion once more.

Nothing in the proofs of this chapter required the Hamiltonian to be independent of time. Step 1 of the proof of the Liouville equation, Eq. (26), only involves the mixed second derivatives of \(\mathcal{H}\), which are equal for any \(\mathcal{H}\). Therefore, the Liouville equation holds and the phase space volume is conserved, but only in the canonical pair \((q,p)\). In terms of the position and velocity, \(v = \dot{q}\), the equations of motion are \(\dot{q} = v\) and \(\dot{v} = - \gamma v - \omega_0^2 q\), so that

\[ \frac{\partial \dot{q}}{\partial q} + \frac{\partial \dot{v}}{\partial v} = - \gamma \neq 0. \]

Repeating the proof of Eq. (32) in these variables gives \(\frac{\mathrm{d}}{\mathrm{d}t} \det{J} = - \gamma \det{J}\), so a volume of initial conditions in the \((q,v)\) plane shrinks as \(e^{-\gamma t}\). Fig. 7 shows the same ensemble of damped oscillators in both planes. In the \((q,v)\) plane, all the oscillators spiral towards rest at the origin and the square shrinks. In the \((q,p)\) plane, it is squeezed in \(q\) but stretched in \(p = m e^{\gamma t} \dot{q}\), and its area stays constant. Phase space volume is conserved in the canonical coordinates, not in any coordinates we find natural.

Fig. 7 An ensemble of damped oscillators (\(m = \omega_0 = 1\), \(\gamma = 0.1\)) starting from the same distribution as in Fig. 5, in the \((q,v)\) plane (top) and in the \((q,p)\) plane (bottom), with the area of the square. Note the different scales of the two panels.

Equilibrium

Now we go back to the connection we established between the Poisson bracket of the phase space distribution function with the Hamiltonian, through Eq. (31). We define thermodynamic equilibrium as the case when the distribution \(f(\textbf{x})\) does not depend explicitly on time. Eq.(31) then reads

\[ \{ f, \mathcal{H} \} = 0. \]

This equation can be solved for \(f(\textbf{x})\) as a function of the Hamiltonian, which we formally denote as

\[ f(\textbf{x}) = \frac{1}{\mathcal{Z}} \mathcal{F}(\mathcal{H}(\textbf{x})). \]

The averages of macroscopic observables will then also be time-independent (you might remember this definition from earlier classes for thermodynamic equilibrium), given by

(35)\[ A = \langle a \rangle_f = \frac{1}{\mathcal{Z}} \int \mathcal{F}(\mathcal{H}(\textbf{x})) \, a(\textbf{x}) \, \mathrm{d} \textbf{x}. \]

The denominator \(\mathcal{Z}\) is called the partition function, which is used to normalize the function \(\mathcal{F}\), defined as

(36)\[ \mathcal{Z} = \int \mathcal{F}(\mathcal{H}(\textbf{x})) \, \mathrm{d} \textbf{x}. \]

To proceed and write specific expressions for \(\mathcal{F}\), we invoke a second postulate, called the hypothesis of equal a priori probabilities. In simple terms, it says that for an ensemble of isolated systems, each characterized by the same number of particles, volume, and total energy, all microstates have equal probability. That is, the phase space distribution is a constant when the energy is between \(E\) and \(E + \delta E\) and zero otherwise.

We saw previously that the classical equations of motion for an isolated system under the influence of conservative forces ensure that the total energy is constant. An ensemble composed of such isolated systems is called the microcanonical ensemble. Since they are characterized by a constant number of particles, volume, and total energy, the ensemble is often denoted by the letters NVE.

According to our second postulate, the equilibrium distribution at the NVE ensemble can be written formally as

\[ \mathcal{F}(\mathcal{H}(\textbf{x})) \propto \delta (\mathcal{H}(\textbf{x}) - E), \]

and expectation values are obtained through

(37)\[ A = \langle a \rangle_E = \frac{ \int \delta (\mathcal{H}(\textbf{x}) - E) \, a(\textbf{x}) \, \mathrm{d} \textbf{x}} {\int \delta (\mathcal{H}(\textbf{x}) - E) \, \mathrm{d} \textbf{x}} . \]

You have derived in previous courses \(\mathcal{F}(\mathcal{H}(\textbf{x}))\) for other ensembles. We will not repeat this derivation here, but just mention that all of them can be derived from our two postulates by considering the entire ensemble as an isolated system. For an ensemble of systems with a fixed number of particles, volume, and temperature, called the canonical ensemble, \(\mathcal{F}\) is the famous Boltzmann distribution,

\[ \mathcal{F}(\mathcal{H}(\textbf{x})) \propto e^{-\beta \mathcal{H} (\textbf{x})}. \]

and the observables are given by

(38)\[ A = \langle a \rangle_{\beta} = \frac{ \int e^{-\beta \mathcal{H} (\textbf{x})} \, a(\textbf{x}) \, \mathrm{d} \textbf{x}} {\int e^{-\beta \mathcal{H} (\textbf{x})} \, \mathrm{d} \textbf{x}} . \]

The canonical ensemble is usually denoted by the letters NVT.

Fig. 8 and Fig. 9 show two examples, an ensemble of harmonic oscillators and an ensemble of pendulums, both in the canonical ensemble. Compare them to Fig. 5 and Fig. 6: the systems move along their trajectories exactly as before, but the distributions do not change in time. Note that we never solved the equations of motion of the pendulum, and we do not need to: \(e^{-\beta \mathcal{H}}\) solves the Liouville equation for any \(\mathcal{H}\).

Fig. 8 An ensemble of harmonic oscillators at equilibrium, \(f \propto e^{-\beta \mathcal{H}}\) (color), with a few systems shown as points and four of them with their trajectories. The systems move, the distribution does not.

Fig. 9 The same for an ensemble of pendulums, over the full circle \(-\pi \le \phi \le \pi\). Systems that go over the top leave the frame on one side and reappear on the other.