thermodynamics - Derivation of the virial expansion

It is easy to derive the ideal gas law form the kinetic theory of gases, but how do you derive the other coefficients of the virial expansion? Is it even possible to get a closed-form formula for the coefficients just from the kinetic theory of gases?

Answer

I am going to show how to get a closed form for the second virial coefficient $B_2(T)$ for a monoatomic gas, by using a potential energy between two atoms at a distance $r$ from each other given by

$$\phi(r) = \begin{cases} +\infty,& r < \sigma \\ -k\epsilon\left(\dfrac{\sigma}{r}\right)^6, & r > \sigma \end{cases}$$

where $k$ is Boltzmann constant, and where $\sigma$ and $\epsilon$ are parameters to be adjusted to fit experimental data: $\sigma$ shall be of the order of the atom radius, whereas $\epsilon$ embodies the strength of the interaction. This is called the Sutherland's potential. The $r^{-6}$ law at large distance is motivated by the fact that atoms are neutral and therefore that their interaction takes place as dipole-dipole (van der Waals forces). The infinitely repulsive law at short distance models the Fermi exclusion principle. You may have seen a law in $r^{-12}$ instead (the Lennard-Jones potential): it tends to fit experimental data slightly better but the computations will already be complicated enough that I prefer to make the short-range side simple at least.

First, let me showcase how well it works, taking Argon as an example: the following graph shows experimental points [GO64][*] in red and the theoretical curve in blue. The experimental error is $\pm 2$ at most in the given units.

Now, you should be impressed and therefore motivated to read the derivation ;-)

First let's start with thermodynamics. Let $N$ be the number of atoms in the gas, occupying a volume $V$, at pressure $P$, and temperature $T$. The virial expansion reads

$$\frac{P}{kT} = \frac{N}{V} + B_2(T)\left(\frac{N}{V}\right)^2 + \cdots \tag{1}\label{virial}$$

We will extract $P$ from

$$P = -\left.\frac{\partial F}{\partial V}\right|_{T,N},$$

where $F$ is the free energy.

Now let's move to statistical physics to compute $F$. I will assume the reader to know basic statistical physics, specifically the central role played by the partition function $Z$, and especially,

$$F=-kT\ln Z.$$

Since $N$ is constant, we are in the case of the canonical ensemble, thus

$$Z = \int \exp\left(-\frac{u_N}{kT}\right)d\tau_N,$$

where $u_N$ is the total energy density for the $N$ atoms. It is the sum of their kinetic energy and of the potential energies of the interaction between pairs of atoms. Let $p_i$ denote the momentum of the $i$-th atom, and $x_i$ its position (both are vectors with 3 component), and $r_{ij}$ the distance between the $i$-the and the $j$-th atom. Then

$$u_N = \sum_{i=1}^N \frac{p_i^2}{2m} + \sum_{i=1}^N \sum_{j=i+1}^N \phi(r_{ij}).$$

As for $d\tau_N$, it is the infinitesimal element of phase space,

$$d\tau_N=\frac{1}{N!}d^3p_1d^3x_1\cdots d^3p_Nd^3x_N.$$

The factorial accounts for the indistinguishability of the atoms but this is not important here as we shall see very soon. Thus, in the expression of $Z$, the integration for each component of momentum ranges from $-\infty$ to $+\infty$, whereas each integration $\int d^3x_i$ shall cover the entire volume $V$.

We see that $u_N$ is the sum of one term depending only on the momenta and of another term depending only on the positions. This allows to factorise the expression of $Z$ as follow,

$$Z=\frac{1}{N!}\int \exp\left(-\frac{1}{kT}\sum_i \frac{p_i^2}{2m}\right)d^3p_1\cdots d^3p_N\int \exp\left(-\sum_{i

For a perfect gas, we would have $\phi=0$, and therefore the partition function would read

$$Z_\text{ideal}=\frac{V^N}{N!}\int \exp\left(-\frac{1}{kT}\sum_i \frac{p_i^2}{2m}\right)d^3p_1\cdots d^3p_N.$$

Therefore

$$F = \underbrace{-kT\ln Z_\text{ideal}}_{F_\text{ideal}} -kT\ln Z_\text{interaction},$$

where

$$Z_\text{interaction}=\frac{1}{V^N}\int \exp\left(-\sum_{i

Of course, we shall recover the perfect gas law from the first term,

$$P=-\frac{\partial F_\text{ideal}}{\partial V} = \frac{NkT}{V}.$$

This is the derivation the OP mentioned in his question: I will assume this result and move on to the second term,

$$P = \frac{NkT}{V} + kT\frac{\partial \ln Z_\text{interaction}}{\partial V},$$

and therefore with the definition of the virial coefficient from eqn (\ref{virial}),

$$B_2(T) = \left(\frac{V}{N}\right)^2 \frac{\partial \ln Z_\text{interaction}}{\partial V}. \tag{2}\label{B2}$$

Noticing that

$$ \exp\left( -\sum_{i

the trick to compute $Z_\text{interaction}$ is to introduce

$$f_{ij} = \exp\left(-\frac{\phi(r_{ij})}{kT}\right) - 1.$$

This quantity is very small except when the two atoms are close to each others because the potential $\phi$ decreases quickly when the distance between the atoms increases. But the virial expansion is an expansion in the density $N/V$: i.e. we assume that $N/V$ is small, and therefore that the atoms are far and apart. Therefore the $f_{ij}$'s are small and we can therefore perform an expansion. So $Z$ now features a term

$$\prod_{i

The third term is small unless at least three atoms are close to each other: it is much smaller than the second term with our low density hypothesis. Higher-order terms are even more negligible. So we will keep only the first two terms. But it should be noted that the third term would contribute toward the 3rd virial coefficients. There is a big theory called cluster expansion, which provides a formal framework to perform this expansion correctly but I won't discuss that.

Thus we have now,

$$Z_\text{interaction} = \frac{1}{V^N}\int\left(1 + \sum_{i

By performing the integration over the positions that $f_{ij}$ does not depend upon, we get

$$Z_\text{interaction} = 1 + \frac{1}{V^2}\sum_{i

Then we take the logarithm of $Z_\text{interaction}$ and proceed with another expansion,

$$\ln(1+w) = w + \mathcal{O}(w^2)$$

in the same spirit as the previous ones,

$$\ln Z_\text{interaction} = \frac{1}{V^2}\sum_{i

But then all these integrals are equal and there are $N(N-1)/2$ of them, which is $N^2/2$ in the large $N$ limit where we necessarily work. Furthermore, $f_{ij}$ depends only on $r_{ij}$, and we can therefore perform one more integration, on $x_i$ e.g., which gives another factor $V$, and then perform the integration on the separation $x_j - x_i$ in spherical coordinates,

$$\ln Z_\text{interaction}=\frac{4\pi N^2}{2V}\int_0^{+\infty} \left(\exp\left(-\frac{\phi(r)}{kT}\right) - 1 \right)r^2 dr,$$

where $4\pi$ comes from the integration over the solid angle.

Thus with eqn (\ref{B2}),

$$B_2(T) = -2\pi\int_0^{+\infty} \left(\exp\left(-\frac{\phi(r)}{kT}\right) - 1 \right)r^2 dr.$$

The computation of this integral is straightforward if boring, in two steps, first the easy one,

$$B_2(T) = \frac{2\pi \sigma^3}{3}\left(1 - 3\int_1^{+\infty}s^2\left(\exp\left(\frac{\epsilon}{T}\frac{1}{s^6}\right)-1\right)ds\right),$$

and then, the guru step,

$$B_2(T) = \frac{2\pi \sigma^3}{3} \Lambda \left(\frac{\epsilon}{T}\right),$$

where

$$\Lambda(w) = \exp w - \sqrt{\pi} w \; \mathrm{erfi}(x).$$

$\mathop{\rm erfi}$ is the so-called imaginary error function, available in Matlab, Mathematical, etc. Then I did a fit of that expression to the experimental data to find $\epsilon$ and $\sigma$. The curve I showed at the beginning corresponds to $\sigma=2.429\ Å$ and $\epsilon=349.42\ \mathrm{K}$.

That's all folks!

[*] Those data are pretty old but I had them lying around in digital form: anyone, please feel free to point me to more recent ones!

[GO64] Donald A. Gyorog and Edward F. Obert. Virial coefficients for argon, methane, nitrogen, and xenon. AIChE Journal, 10(5):621–625, 1964.