Virial Expansion for Hard Spheres

In this post I derive the virial expansion for a classical gas of hard spheres. The nice part is watching the excluded volume correction come out of the algebra on its own, and at the end I look at the higher-order coefficients (which lead to a surprisingly good closed-form equation of state).

What the virial expansion actually computes#

The ideal gas law $PV = Nk_BT$ works well, but it makes two simplifying assumptions that break down for real gases:

Particles are point-like (zero volume)
Particles don’t interact (no forces between them)

Real atoms and molecules violate both: they have finite size and exert forces on each other. The virial expansion systematically corrects for these effects by expanding in powers of the density:

\frac{PV}{Nk_BT} = 1 + B_2\rho + B_3\rho^2 + \cdots

where $\rho = N/V$ is the number density. The virial coefficients $B_2, B_3, \ldots$ encode how interactions and finite size change the pressure. At low density ( $\rho \to 0$ ), we recover the ideal gas. As density increases, higher-order terms become important.

The hard sphere model captures the simplest finite-size effect: particles are impenetrable balls of diameter $\sigma$ . The potential is:

U(r) = \begin{cases} \infty & \text{if } r < \sigma \\ 0 & \text{if } r > \sigma \end{cases}

This is more than a toy model (it is the starting point for real gases, liquids and even the glass transition). So the question is: how does this simple exclusion change the equation of state?

The partition function and why it’s hard to compute#

For $N$ particles with Hamiltonian:

H = \sum_{i=1}^N \frac{p_i^2}{2m} + \sum_{i>j} U(r_{ij})

the classical partition function separates into momentum and position integrals:

Z(N, V, T) = \frac{1}{N!\lambda^{3N}} \int \prod_i d^3r_i\, e^{-\beta\sum_{j<k} U(r_{jk})}

where $\lambda = \sqrt{2\pi\hbar^2/(mk_BT)}$ is the thermal wavelength. The momentum integral is easy (it gives the ideal gas result), but the position integral is the challenge:

Q(N, V, T) = \int \prod_i d^3r_i\, e^{-\beta\sum_{j<k} U(r_{jk})}

For interacting particles, the Boltzmann factor $e^{-\beta U}$ doesn’t factorize, so we can’t separate this integral into $N$ independent pieces. This is why we need the Mayer expansion.

The Mayer f-function: turning a hard problem into a series#

The trick is to define:

f(r) = e^{-\beta U(r)} - 1

For hard spheres, this is very simple:

f(r) = \begin{cases} -1 & \text{if } r < \sigma \\ 0 & \text{if } r > \sigma \end{cases}

With this, the partition function becomes:

Q = \int \prod_i d^3r_i\, \prod_{j>k} (1 + f_{jk})

where $f_{jk} = f(r_{jk})$ . Now we can expand the product:

\prod_{j>k} (1 + f_{jk}) = 1 + \sum_{j>k} f_{jk} + \sum_{j>k, l>m} f_{jk}f_{lm} + \cdots

Each term represents a distinct cluster of interacting particles. The first term ( $1$ ) gives the ideal gas. The linear terms ( $\sum f_{jk}$ ) represent two-particle interactions. The quadratic terms represent three- and four-particle clusters, and so on.

Physically, the Mayer $f$ -function measures how far the Boltzmann factor of a pair deviates from the non-interacting value of 1. $f = 0$ means no interaction (particles independent). $f = -1$ means complete exclusion (particles can’t overlap).

Keeping only the first correction: why this is enough at low density#

At low density, most particles are far apart and don’t interact. We only need to keep terms where a small number of particles are close together. Keeping only the linear term in $f$ :

Q \approx \int \prod_i d^3r_i\, \left(1 + \sum_{j>k} f_{jk}\right)

The first term gives $V^N$ (the ideal gas contribution). The second term requires evaluating:

\int \prod_i d^3r_i\, f_{12} = V^{N-2} \int d^3r_1\,d^3r_2\, f(r_{12})

All $N(N-1)/2$ pairs contribute equally, so we just need to compute one integral and multiply by the number of pairs.

Now change to center-of-mass and relative coordinates:

\vec{R} = \frac{1}{2}(\vec{r}_1 + \vec{r}_2), \qquad \vec{r} = \vec{r}_1 - \vec{r}_2

The Jacobian is $1$ , and the integral separates:

\int d^3r_1\,d^3r_2\, f(r_{12}) = \left(\int d^3R\right) \left(\int d^3r\, f(r)\right) = V \int d^3r\, f(r)

The center-of-mass integral just gives the volume $V$ . The relative coordinate integral is what contains the physics.

The second virial coefficient: what it actually means#

Carrying through the algebra, the equation of state becomes:

\frac{PV}{Nk_BT} = 1 - \frac{N}{2V}\int d^3r\, f(r)

For hard spheres:

\int d^3r\, f(r) = -4\pi \int_0^\sigma r^2\, dr = -\frac{4\pi\sigma^3}{3}

This integral is the volume excluded by one sphere (of radius $\sigma$ ), but with a negative sign because $f = -1$ in the excluded region.

Plugging this in:

\frac{PV}{Nk_BT} = 1 + \frac{2\pi\sigma^3}{3}\rho = 1 + B_2\rho

where the second virial coefficient is:

\boxed{B_2 = \frac{2\pi\sigma^3}{3} = 4 \times \frac{4\pi}{3}\left(\frac{\sigma}{2}\right)^3}

Geometrically, $B_2$ is four times the volume of one sphere. Where does the 4 come from? Each particle excludes a sphere of radius $\sigma$ around itself (volume $\frac{4\pi}{3}\sigma^3$ ), and we have to divide by 2 to avoid double-counting pairs.

Excluded volume: the physical picture#

Rewriting the result as:

\frac{P}{k_BT} = \frac{N}{V - V_{\text{ex}}}

where $V_{\text{ex}} = N \times \frac{2\pi\sigma^3}{3}$ is the excluded volume, we get a very clean statement:

The hard sphere gas behaves like an ideal gas with reduced volume.

Two spheres of diameter $\sigma$ can’t have their centers closer than $\sigma$ , so each particle forbids the center of any other particle from a sphere of radius $\sigma$ around itself. That region has volume $\frac{4\pi}{3}\sigma^3$ , which is 8 times the particle’s own volume $\frac{\pi\sigma^3}{6}$ . Each such exclusion is shared by a pair, so each particle is charged half of it: $\frac{2\pi}{3}\sigma^3$ per particle, i.e. 4 particle volumes. For $N$ particles this gives $V_{\text{ex}} = \frac{2}{3}\pi N\sigma^3$ .

So the pressure is higher than the ideal gas at the same density and temperature:

P_{\text{hard sphere}} > P_{\text{ideal}}

This makes sense: collisions between finite-sized particles transfer momentum more effectively than non-interacting point particles. The particles “bounce off” each other, increasing the pressure on the container walls.

Higher-order corrections#

The virial expansion doesn’t stop at $B_2$ . The next term involves three-particle clusters and gives $B_3$ . For hard spheres:

B_2 \propto \sigma^3, \quad B_3 \propto \sigma^6, \quad B_4 \propto \sigma^9, \quad \ldots

Each higher coefficient involves more complicated cluster integrals. A convenient variable is the packing fraction $x = \frac{\pi}{6}\rho\sigma^3$ , the fraction of the volume actually occupied by the spheres (equivalently $x = V_{\text{ex}}/4V$ ):

\frac{PV}{Nk_BT} = 1 + \sum_{n=1}^\infty a_n x^n

where the coefficients $a_n$ are pure numbers determined by geometry. The exact values are:

a_1 = 4,\quad a_2 = 10,\quad a_3 \approx 18.36,\quad a_4 \approx 28.22,\quad a_5 \approx 39.8,\quad \ldots

Only $a_1$ (from $B_2$ ) and $a_2$ (from $B_3$ ) are integers. $a_3$ is known in closed form, but it’s an ugly expression involving $\arccos(1/3)$ . The higher ones come from increasingly nasty cluster integrals, evaluated numerically. They lie very close to the integers $n^2 + 3n$ , i.e. $4, 10, 18, 28, 40, \ldots$ . Summing that integer sequence exactly gives the famous Carnahan–Starling equation of state:

\boxed{\frac{PV}{Nk_BT} \approx 1 + \frac{2(-2 + x)x}{(-1 + x)^3} = \frac{1 + x + x^2 - x^3}{(1 - x)^3}}

This is an approximation, not an exact resummation, but an excellent one: it matches simulations of the hard-sphere fluid to within about a percent all the way up to the freezing point, $x \approx 0.49$ . The pole at $x = 1$ is unphysical. Hard spheres can’t be packed beyond $x = \pi/\sqrt{18} \approx 0.74$ (close packing), and the fluid branch ends well before that.

When does this matter? The expansion parameter is $x \sim \rho\sigma^3$ . For gases at STP, $\rho\sigma^3 \sim 10^{-3}$ , so $B_2$ is a small correction. For liquids, $\rho\sigma^3 \sim 1$ , and many terms are needed.

Where this leaves us#

So the virial expansion is a systematic way to go from point particles to a gas of particles with finite size, and already the first correction is the excluded volume: the particles have less room to move than in an ideal gas, so they hit the walls more often and the pressure goes up. For any other intermolecular potential you do the same thing, just with a different $f(r)$ in the cluster integrals.