How is the van der Waals equation a mean-field model? The van der Waals equation for real gases is
$$\left(  P+a\frac{n^{2}}{V^{2}}\right)\left(  V-nb\right)=nRT$$
where $a$ and $b$ are constants, $n$ is the number of moles, $R$ is the gas constant, and $P$, $T$, and $V$ are the pressure, temperature, and volume, respectively.
The first term in parentheses accounts for the attractive interactions, the second term accounts for the repulsive interactions. I have heard that this is a "mean-field" model - how is this the case? I assume that is because the attractive interactions do not seem to depend on a particular atomic nature, though the repulsive interactions certainly do. What is the utility of the mean-field assumption, and are there non-mean-field models of liquid/vapor equilibria that contrast with the van der Waals model?
 A: Statistical mechanics
Van der Waals equation can be derived from virial/cluster expansion in terms of the average/mean density of the liquid/gas:
$$
\rho(\mathbf{r})=\left\langle \sum_{i=1}^N\delta(\mathbf{r}-\mathbf{r}_i)\right\rangle,
$$
where $\mathbf{r}_i$ is the position of the i-th atom/molecule.
Thermodynamics
Another way to look at it is that pressure is the derivative of the free energy in respect to volume, so one can obtain the expression for the free energy by integrating the pressure expressed from van der Waals equation in respect to volume. One can then analyze liquid-gas transition in terms of Landau theory of phase transitions, with density serving as the order parameter.
See, e.g., section 8.5 in Lecture notes on Statistical Physics by Lydéric Bocquet
See also Wikipedia for the expression for the Helmholz free energy in more conventional terms.
Simple mean-field derivation of Van der Waals equation
Simple mean-field derivation of VdW equation can be performed along the lines of the derivation given in Wikipedia. I sketch it here in even more simplified/clarified form.
The energy of a gas of N interacting particles is
$$
H(\mathbf{x}_1, \mathbf{p}_1;...;\mathbf{x}_N, \mathbf{p}_N)=
\sum_{i=1}^N\frac{\mathbf{p}_i^2}{2m}+\frac{1}{2}\sum_{i=1}^N\sum_{j=1,j\neq i}^NU(|\mathbf{x}_i-\mathbf{x}_j|),
$$
where the potential is attractive, i.e., $U(|\mathbf{x}_i-\mathbf{x}_j|)<0$ (unless the molecules/atoms are very close to each other)
The partition function is
$$
Z=\int \prod_{i=1}^Nd^3\mathbf{x}_id^3\mathbf{p}_ie^{-\beta H(\mathbf{x}_1, \mathbf{p}_1;...;\mathbf{x}_N, \mathbf{p}_N)}=
\int d^3\mathbf{p}_1...d^3\mathbf{p}_Ne^{-\beta \sum_{i=1}^N\frac{\mathbf{p}_i^2}{2m}}Z_c=\\\left(2\pi mk_B T\right)^{\frac{3N}{2}}Z_c
$$
where the configuration integral is
$$
Z_c=\int d^3\mathbf{x}_1...d^3\mathbf{x}_Ne^{-\frac{\beta}{2}\sum_{i=1}^N\sum_{j=1,j\neq i}^NU(|\mathbf{x}_i-\mathbf{x}_j|)},
$$
and $\beta=1/(k_BT)$.
We now perform the mean-field approximation replacing the pairwise interactions between particles by an interaction of each particle with an averaged field of all other particles:
$$
\frac{1}{2}\sum_{i=1}^N\sum_{j=1,j\neq i}^NU(|\mathbf{x}_i-\mathbf{x}_j|)\approx
\frac{1}{2}\sum_{i=1}^N \frac{1}{V}\int_V d^3\mathbf{x_j}U(|\mathbf{x}_i-\mathbf{x}_j|)\approx\\
\frac{1}{2}\sum_{i=1}^N \frac{N}{V}\int_V d^3\mathbf{y}U(|\mathbf{y}|)=-\frac{a'N^2}{V},
$$
where the minus sign is due to the fact that the potential is attractive, so that $a'>0$.
The configuration integral is then
$$
Z_c=\int d^3\mathbf{x}_1...d^3\mathbf{x}_N e^{\beta\frac{a'N^2}{V}}=
(V-b'N)^N e^{\beta\frac{a'N^2}{V}},
$$
where we took account for the fact that the particles cannot be too close to each other, and hence the volume integral for each particle is reduced by the factor accounting for the volume of all other particles.
We now can write down the free energy of the gas as
$$
F=-\frac{1}{\beta}\log Z=-\frac{a'N^2}{V}-\frac{N}{\beta}\log(V-b'N)-\frac{3N}{2\beta}\log\left(2\pi mk_BT\right)
$$
The gas pressure is given by (since $dF=-SdT-pdV$)
$$
P=-\left(\frac{\partial F}{\partial V}\right)_T=
\frac{Nk_BT}{V-Nb'}-\frac{a'N^2}{V^2}
$$
Converting to molar fractions $n=N/N_A$, $b=N_Ab'$, $a=a'N_A^2$, we obtain the usual form of the VdW equation
$$
P=\frac{nRT}{V-nb}-\frac{an^2}{V^2}\Leftrightarrow
\left(P+\frac{an^2}{V^2}\right)(V-nb)=nRT.
$$
