# Mean field theories for weakly interacting systems

If a system is composed of weakly interacting particles we can use perturbative methods to compute the partition function of such systems. For example, consider a fluid composed of ${\displaystyle N}$ particles in a region of volume ${\displaystyle V}$, interacting through a generic two-body potential that depends only on the relative distance between the particles:

${\displaystyle U(\lbrace {\vec {r}}\rbrace )={\frac {1}{2}}\sum _{i\neq j=1}^{N}\varphi (|{\vec {r}}_{i}-{\vec {r}}_{j}|)}$
Its Hamiltonian will be:
${\displaystyle {\mathcal {H}}=\sum _{i=1}^{N}{\frac {{\vec {p}}_{i}{}^{2}}{2m}}+{\frac {1}{2}}\sum _{i\neq j=1}^{N}\varphi (|{\vec {r}}_{i}-{\vec {r}}_{j}|)}$
and its partition function:
${\displaystyle Z_{N}={\frac {1}{h^{3N}N!}}\int \prod _{i=1}^{N}d{\vec {q}}_{i}d{\vec {p}}_{i}e^{-\beta {\mathcal {H}}}={\frac {1}{N!\Lambda ^{3N}}}Q_{N}}$
where ${\displaystyle \Lambda }$ is the thermal wavelength and ${\displaystyle Q_{N}}$ is the configurational sum:
${\displaystyle Q_{N}=\int \prod _{i}d{\vec {r}}_{i}e^{-\beta U(\lbrace {\vec {r}}_{i}\rbrace )}}$
Of course for ideal gases ${\displaystyle U=0}$ and so ${\displaystyle Q_{N}=V^{N}}$, and the dependence on the temperature is included only in ${\displaystyle \Lambda }$; if we consider also the interaction terms we must insert a correction ${\displaystyle \chi }$ in the configurational contribution to the partition function[1]:
${\displaystyle Q_{N}=V^{N}\cdot \chi (N,V,T)}$
which (depending on the possible presence of attractive terms in the interaction potential ${\displaystyle \varphi }$) can in general be also a function of the temperature ${\displaystyle T}$; furthermore the correction depends strongly on the gas density: if it is low the particles will not "perceive" the presence of the other ones and the ideal gas approximation is a good one, while for high densities the particles will be closer to each other and corrections to ${\displaystyle Q_{N}}$ are necessary. Let us note that inserting the correction ${\displaystyle \chi }$ the free energy of the system will be:
${\displaystyle F=F_{\text{ideal}}-k_{B}T\ln \chi }$

## Virial and cluster expansion

The virial expansion is a systematic approach that can be used to incorporate the corrections due to the interactions between the particles, and as we will shortly see it can be obtained from a much general method, the cluster expansion. The virial expansion consists in expanding the thermodynamic quantities of a system in powers of the density ${\displaystyle \rho }$; for example, the virial expansion for the pressure of a gas is:

${\displaystyle {\frac {P}{k_{B}T}}=\rho +B_{2}\rho ^{2}+B_{3}\rho ^{3}+\cdots }$
where ${\displaystyle B_{n}}$ are called virial coefficients, and in general they can depend on the temperature. The virial expansion is very useful because the coefficients ${\displaystyle B_{n}}$ can be experimentally measured (for example, in the case of the pressure they can be determined by properly fitting the isotherms of a system), and as we will see they can be related to microscopic properties of the interparticle interaction. Let us see for example the virial expansion of the Van der Waals equation. From Van der Waals equation we have:
${\displaystyle {\frac {P}{k_{B}T}}={\frac {N}{V-Nb}}-{\frac {aN^{2}}{k_{B}TV^{2}}}}$
which can be rewritten as:
${\displaystyle {\frac {P}{k_{B}T}}={\frac {N}{V}}\left(1-b{\frac {N}{V}}\right)^{-1}-{\frac {a}{k_{B}T}}\left({\frac {N}{V}}\right)^{2}=\rho (1-b\rho )^{-1}-{\frac {a}{k_{B}T}}\rho ^{2}}$
where of course ${\displaystyle \rho =N/V}$ is the density of the system. Therefore, expanding the first term in a Taylor series:
${\displaystyle {\frac {P}{k_{B}T}}=\rho +\left(b-{\frac {a}{k_{B}T}}\right)\rho ^{2}+b^{2}\rho ^{3}+b^{3}\rho ^{4}+\cdots }$
We can thus immediately identify the first virial coefficient:
${\displaystyle B_{2}(T)=b-{\frac {a}{k_{B}T}}}$
The temperature ${\displaystyle T_{B}}$ at which ${\displaystyle B_{2}}$ vanishes is in general called Boyle temperature; in this case ${\textstyle T_{B}={\frac {a}{bk_{B}}}}$.

Now, let us see how the cluster expansion works and how we can obtain the virial expansion from it. Of course, we start from the general configurational partition function:

${\displaystyle Q_{N}=\int \prod _{i}d{\vec {r}}_{i}e^{-\beta U(\lbrace {\vec {r}}_{i}\rbrace )}}$
The idea is to find a "small quantity" in terms of which we can expand ${\displaystyle Q_{N}}$; this quantity is the so called Mayer function:
${\displaystyle f({\vec {r}})=e^{-\beta \varphi ({\vec {r}})}-1}$
In fact, when the gas is ideal ${\displaystyle f({\vec {r}})=0}$, and if the particles interact weakly ${\displaystyle \varphi }$ is small, and so is ${\displaystyle f({\vec {r}})}$. In particular, this expansion will work well for low densities (namely ${\displaystyle |{\vec {r}}_{i}-{\vec {r}}_{j}|}$ is large and so ${\displaystyle \varphi (|{\vec {r}}_{i}-{\vec {r}}_{j}|)\to 0}$) or high temperatures (namely ${\displaystyle \beta \to 0}$): in both cases, in fact, ${\displaystyle e^{-\beta \varphi (|{\vec {r}}_{i}-{\vec {r}}_{j}|)}\to 1}$ and ${\displaystyle f(|{\vec {r}}_{i}-{\vec {r}}_{j}|)\to 0}$. Using the short notations ${\displaystyle \varphi _{ij}=\varphi (|{\vec {r}}_{i}-{\vec {r}}_{j}|)}$ and ${\displaystyle f_{ij}=f(|{\vec {r}}_{i}-{\vec {r}}_{j}|)}$ we have:
{\displaystyle {\begin{aligned}e^{-\beta {\frac {1}{2}}\sum _{i\neq j=1}^{N}\varphi _{ij}}=\prod _{i\neq j}e^{-{\frac {\beta }{2}}\varphi _{ij}}=\prod _{i\neq j}{\sqrt {1+f_{ij}}}=\\={\sqrt {1+f_{12}}}{\sqrt {1+f_{13}}}\cdots {\sqrt {1+f_{1N}}}{\sqrt {1+f_{21}}}{\sqrt {1+f_{23}}}\cdots \\\cdots {\sqrt {1+f_{N1}}}{\sqrt {1+f_{N2}}}\cdots {\sqrt {1+f_{N,N-1}}}\end{aligned}}}
Now, since ${\displaystyle f_{ij}=f_{ji}}$ by definition, in the end we will have:
${\displaystyle e^{-\beta {\frac {1}{2}}\sum _{i\neq j=1}^{N}\varphi _{ij}}=(1+f_{12})(1+f_{13})\cdots =1+{\frac {1}{2}}\sum _{i\neq j=1}^{N}f_{ij}+\cdots }$
where the missing terms contain the product of two or more Mayer functions. Therefore, the configurational contribution to the partition function will be:
{\displaystyle {\begin{aligned}Q_{N}=\int \prod _{k=1}^{N}d{\vec {r}}_{k}\left(1+{\frac {1}{2}}\sum _{i\neq j=1}^{N}f_{ij}+\cdots \right)=V^{N}+\int \prod _{k=1}^{N}d{\vec {r}}_{k}{\frac {1}{2}}\sum _{i\neq j=1}^{N}f_{ij}+\cdots =\\=V^{N}+V^{N-2}{\frac {1}{2}}\sum _{i\neq j=1}^{N}\int d{\vec {r}}_{i}d{\vec {r}}_{j}f_{ij}+\cdots \end{aligned}}}
where in the last step we have extracted the contributions of the integrals with ${\displaystyle k\neq i,j}$. Now, the remaining integral can be easily computed with the definition of the new variable ${\displaystyle {\vec {r}}={\vec {r}}_{i}-{\vec {r}}_{j}}$:
${\displaystyle \int d{\vec {r}}_{i}d{\vec {r}}_{j}f_{ij}=\int d{\vec {r}}_{i}d{\vec {r}}_{j}f(|{\vec {r}}_{i}-{\vec {r}}_{j}|)=\int d{\vec {r}}_{i}\int d{\vec {r}}f(|{\vec {r}}|)=V\int d{\vec {r}}f(|{\vec {r}}|):=-2VB_{2}}$
where in the last step we have defined the first virial coefficient. Of course, we should have set ${\displaystyle \int d{\vec {r}}f(|{\vec {r}}|)}$ equal to a generic coefficient, but in the end we should have found exactly that this coefficient is ${\displaystyle -2B_{2}}$: we have done this for the sake of simplicity. Therefore:
${\displaystyle B_{2}=-{\frac {1}{2}}\int f(|{\vec {r}}|)d{\vec {r}}}$
From this we see precisely how the virial coefficient, which as we have already stated can be experimentally measured, is related to the microscopic properties of the interaction between the particles, represented by the Mayer function ${\displaystyle f}$. It can also be shown that all the virial coefficients can be expressed in terms of integrals of products of Mayer functions; for example the second virial coefficient is:
${\displaystyle B_{3}=-{\frac {1}{3}}\int d{\vec {r}}d{\vec {s}}f(|{\vec {r}}|)f(|{\vec {s}}|)f(|{\vec {r}}-{\vec {s}}|)}$
Higher order coefficients involve the computation of increasingly difficult integrals, which can however be visualized in terms of graphs.

What we have seen now is how the cluster expansion works in general. Let us now apply it in order to find the virial expansion for real gases. From what we have found, the configurational partition function of the system becomes:

${\displaystyle Q_{N}=V^{N}-V^{N-1}B_{2}\sum _{i\neq j=1}^{N}1+\cdots }$
The remaining sum is equal to ${\displaystyle N(N-1)}$: in fact, for any of the ${\displaystyle N}$ values that ${\displaystyle i}$ can assume, ${\displaystyle j}$ can have ${\displaystyle N-1}$ values. Therefore:
${\displaystyle Q_{N}=V^{N}-V^{N-1}N(N-1)B_{2}+\cdots }$
and, considering that ${\displaystyle N-1\approx N}$ for large ${\displaystyle N}$, the complete partition function of the system will be:
${\displaystyle Z_{N}={\frac {V^{N}}{N!\Lambda ^{3N}}}\left(1-{\frac {N^{2}}{V}}B_{2}+\cdots \right)}$
We recognise in this expression that ${\displaystyle (1+B_{2}N^{2}/V+\cdots )}$ is the correction ${\displaystyle \chi }$ to the ideal gas partition function that we have mentioned earlier; therefore, the free energy of the system will be:
${\displaystyle F=F_{\text{ideal}}-k_{B}T\ln \left(1-{\frac {N^{2}}{V}}B_{2}+\cdots \right)}$
and its pressure:
${\displaystyle P=-{\frac {\partial F}{\partial V}}_{|N,T}={\frac {Nk_{B}T}{V}}+k_{B}T{\frac {{\frac {N^{2}}{V^{2}}}B_{2}}{1-{\frac {N^{2}}{V}}B_{2}}}+\cdots ={\frac {Nk_{B}T}{V}}\left(1+{\frac {{\frac {N}{V}}B_{2}}{1-{\frac {N^{2}}{V}}B_{2}}}+\cdots \right)}$
Now, neglecting the terms involving ${\displaystyle B_{n}}$ with ${\displaystyle n\geq 3}$, and expanding in terms of ${\displaystyle B_{2}}$ (in fact ${\displaystyle B_{2}\sim f}$, which is small)[2]:
${\displaystyle P={\frac {Nk_{B}T}{V}}\left(1+{\frac {N}{V}}B_{2}+\cdots \right)}$
This expansion contains only low-order terms in the density ${\displaystyle N/V}$, so strictly speaking it is valid only for low densities. However, we can use a "trick" in order to extend its range; in fact, remembering that the McLaurin expansion of ${\displaystyle (1-x)^{-1}}$ is ${\displaystyle 1+x+\cdots }$, from the last equation we can write:
${\displaystyle {\frac {PV}{Nk_{B}T}}=1+B_{2}\rho +\cdots \approx {\frac {1}{1-B_{2}\rho }}}$
and now re-expand ${\displaystyle (1-B_{2}\rho )^{-1}}$, so that we can express all the virial coefficients in terms of the first one:
{\displaystyle {\begin{aligned}{\frac {1}{1-B_{2}\rho }}=1+B_{2}\rho +(B_{2})^{2}\rho ^{2}+(B_{2})^{3}\rho ^{3}+\cdots \quad \Rightarrow \\\Rightarrow \quad {\frac {PV}{Nk_{B}T}}={\frac {P}{\rho k_{B}T}}=1+B_{2}\rho +(B_{2})^{2}\rho ^{2}+(B_{2})^{3}\rho ^{3}+\cdots \quad \Rightarrow \\\Rightarrow \quad {\frac {P}{k_{B}T}}=\rho +B_{2}\rho ^{2}+(B_{2})^{2}\rho ^{3}+(B_{2})^{3}\rho ^{4}+\cdots \end{aligned}}}
So in the end:
${\displaystyle B_{3}\sim (B_{2})^{2}\qquad B_{4}\sim (B_{2})^{3}\qquad \cdots \qquad B_{n}\sim (B_{2})^{n-1}}$

## Computation of virial coefficients for some interaction potentials

Let us now see this method in action by explicitly computing some coefficients ${\displaystyle B_{2}}$ for particular interaction potentials.

### Hard sphere potential

As a first trial, we use a hard sphere potential similar to the one we have seen for the derivation of the Van der Waals equation:

${\displaystyle \varphi ({\vec {r}})={\begin{cases}+\infty &|{\vec {r}}|
(the difference with what we have seen in Van der Waals equation is that now the potential is purely repulsive, and has no attractive component). In this case:
${\displaystyle f({\vec {r}})=e^{-\beta \varphi ({\vec {r}})}-1={\begin{cases}-1&|{\vec {r}}|
Therefore, from the definition of ${\displaystyle B_{2}}$ and shifting to spherical coordinates:
${\displaystyle B_{2}=-{\frac {1}{2}}\int f(|{\vec {r}}|)d{\vec {r}}=-{\frac {1}{2}}\int _{0}^{+\infty }4\pi r^{2}\left(e^{-\beta \varphi (r)}-1\right)dr=2\pi \int _{0}^{r_{0}}r^{2}dr={\frac {2}{3}}\pi r_{0}^{3}}$
In this case ${\displaystyle B_{2}}$ does not depend on the temperature, and in the end we have:
${\displaystyle PV=Nk_{B}T\left(1+{\frac {2}{3}}\pi r_{0}^{3}{\frac {N}{V}}\right)}$

### Square well potential

We now use a slight refinement of the previous potential:

${\displaystyle \varphi ({\vec {r}})={\begin{cases}+\infty &|{\vec {r}}|r_{0}+\delta \end{cases}}}$
This can be seen as a hard sphere potential where the spheres have an attractive shell of thickness ${\displaystyle \delta }$. We thus have:
${\displaystyle f({\vec {r}})={\begin{cases}-1&|{\vec {r}}|r_{0}+\delta \end{cases}}}$
so that:
{\displaystyle {\begin{aligned}B_{2}=-{\frac {1}{2}}\int f(|{\vec {r}}|)d{\vec {r}}=-{\frac {1}{2}}\int 4\pi r^{2}f(r)dr=\\=-2\pi \left[\int _{0}^{r_{0}}(-r^{2})dr+\int _{r_{0}}^{r_{0}+\delta }\left(e^{\beta \varepsilon }-1\right)r^{2}dr\right]=\\{}\\=-2\pi \left\lbrace -{\frac {r_{0}^{3}}{3}}+{\frac {e^{\beta \varepsilon }-1}{3}}\left[(r_{0}+\delta )^{3}-r_{0}^{3}\right]\right\rbrace =B_{2}^{\text{h.s.}}-{\frac {2}{3}}\pi \left(e^{\beta \varepsilon }-1\right)\left[(r_{0}+\delta )^{3}-r_{0}^{3}\right]\end{aligned}}}
where ${\displaystyle B_{2}^{\text{h.s.}}}$ is the first virial coefficient of the hard sphere potential we have previously seen. Now, if the temperature is sufficiently high, namely ${\displaystyle \beta \varepsilon \ll 1}$, we can approximate ${\displaystyle e^{\beta \varepsilon }-1\approx \beta \varepsilon }$, so that:
${\displaystyle B_{2}=B_{2}^{\text{h.s.}}-{\frac {2}{3}}\pi \beta \varepsilon r_{0}^{3}\left[\left(1+{\frac {\delta }{r_{0}}}\right)^{3}-1\right]}$
For the sake of simplicity, defining:
${\displaystyle \lambda :=\left(1+{\frac {\delta }{r_{0}}}\right)^{3}-1}$
we will have, in the end:
${\displaystyle {\frac {PV}{Nk_{B}T}}=1+B_{2}\rho =1+\left(B_{2}^{\text{h.s.}}-{\frac {2}{3}}{\frac {\pi \varepsilon }{k_{B}T}}r_{0}^{3}\lambda \right)\rho }$
so in this case ${\displaystyle B_{2}}$ actually depends on the temperature.

### Lennard-Jones potential

This potential is a quite realistic representation of the interatomic interactions. It is defined as:

${\displaystyle \varphi ({\vec {r}})=4\varepsilon \left[\left({\frac {r_{0}}{|{\vec {r}}|}}\right)^{12}-\left({\frac {r_{0}}{|{\vec {r}}|}}\right)^{6}\right]}$
which contains a long-range attractive term (the one proportional to ${\displaystyle 1/|{\vec {r}}|^{6}}$, which can be justified in terms of electric dipole fluctuations) and a short-range repulsive one (proportional to ${\displaystyle 1/|{\vec {r}}|^{12}}$, which comes from the overlap of the electron orbitals). With this interaction potential, the first virial coefficient is:
${\displaystyle B_{2}=-2\pi \int _{0}^{+\infty }r^{2}\left\lbrace e^{-{\frac {4\varepsilon }{k_{B}T}}\left[\left({\frac {r_{0}}{|{\vec {r}}|}}\right)^{12}-\left({\frac {r_{0}}{|{\vec {r}}|}}\right)^{6}\right]}-1\right\rbrace dr}$
which is not analytically computable. However, it can be simplified defining the variables ${\displaystyle x=r/r_{0}}$ and ${\displaystyle T^{*}=k_{B}T/\varepsilon }$ so that:
${\displaystyle B_{2}={\frac {2}{3}}\pi r_{0}^{3}{\frac {4}{T^{*}}}\int _{0}^{\infty }x^{2}\left({\frac {12}{x^{12}}}-{\frac {6}{x^{6}}}\right)e^{-{\frac {4}{T^{*}}}\left({\frac {1}{x^{12}}}-{\frac {1}{x^{6}}}\right)}dx}$
Now, we can expand the exponential and integrate term by term; this gives an expression of ${\displaystyle B_{2}}$ as a power series of ${\displaystyle 1/T^{*}}$:
${\displaystyle B_{2}=-2A\sum _{n=0}^{\infty }{\frac {1}{4n!}}\Gamma \left({\frac {2n-1}{4}}\right)\left({\frac {1}{T^{*}}}\right)^{(2n+1)/4}}$
where ${\displaystyle \Gamma }$ is the Euler gamma function and ${\displaystyle A}$ is a constant. Note that the attractive part of the Lennard-Jones potential has introduced in ${\displaystyle B_{2}}$ a dependence on the temperature.

1. We can always insert a multiplicative correction. We could have also written it as an additive correction, but the core of the subject doesn't change.
2. Referring to what we have stated previously, if we set ${\displaystyle \int d{\vec {r}}f(|{\vec {r}}|)={\mathfrak {F}}}$ with ${\displaystyle {\mathfrak {F}}}$ a generic constant, then we should have found:
${\displaystyle Q_{N}=V^{N}\left(1+{\frac {N(N-1)}{V}}\cdot {\frac {\mathfrak {F}}{2}}+\cdots \right)}$
and proceeding like we have done now, in the end:
${\displaystyle P={\frac {Nk_{B}T}{V}}\left(1-{\frac {N}{V}}\cdot {\frac {\mathfrak {F}}{2}}+\cdots \right)}$
and so we see that indeed ${\displaystyle {\mathfrak {F}}=-2B_{2}}$.