**3. The structure of liquids**

### **3.1 Scattered radiation in liquids**

The pair correlation function *g*(*r*) can be deduced from the experimental measurement of the structure factor *S*(*q*) by X-ray, neutron or electron diffractions. In condensed matter, the scatterers are essentially individual atoms, and diffraction experiments can only measure the structure of monatomic liquids such as rare gases and metals. By contrast, they provide no information on the structure of molecular liquids, unless they are composed of spherical molecules or monatomic ions, like in some molten salts.

Furthermore, each type of radiation-matter interaction has its own peculiarities. While the electrons are diffracted by all the charges in the atoms (electrons and nuclei), neutrons are diffracted by nuclei and X-rays are diffracted by the electrons localized on stable electron shells. The electron diffraction is practically used for fluids of low density, whereas the beams of neutrons and X-rays are used to study the structure of liquids, with their advantages and disadvantages. For example, the radius of the nuclei being 10, 000 times smaller than that of atoms, it is not surprising that the structure factors obtained with neutrons are not completely identical to those obtained with X-rays.

To achieve an experience of X-ray diffraction, we must irradiate the liquid sample with a monochromatic beam of X-rays having a wavelength in the range of the interatomic distance (*<sup>λ</sup>* <sup>∼</sup> 0, 1 nm). At this radiation corresponds a photon energy (*h<sup>ν</sup>* <sup>=</sup> *hc <sup>λ</sup>* <sup>∼</sup> <sup>10</sup><sup>4</sup> eV), much larger than the mean energy of atoms that is of the order of few *kBT*, namely about 10−<sup>1</sup> eV. The large difference of the masses and energies between a photon and an atom makes that the photon-atom collision is elastic (constant energy) and that the liquid is transparent to the radiation. Naturally, the dimensions of the sample must be sufficiently large compared to the wavelength *λ* of the radiation, so that there are no side effects due to the walls of the enclosure - but not too much though for avoiding excessive absorption of the radiation. This would be particularly troublesome if the X-rays had to pass across metallic elements with large atomic numbers.

The incident radiation is characterized by its wavelength *λ* and intensity *I*0, and the diffraction patterns depend on the structural properties of the liquids and on the diffusion properties of atoms. In neutron scattering, the atoms are characterized by the scattering cross section *σ* = <sup>4</sup>*πb*2, where *<sup>b</sup>* is a parameter approximately equal to the radius of the core (<sup>∼</sup> <sup>10</sup>−<sup>14</sup> m). Note that the parameter *b* does not depend on the direction of observation but may vary slightly, even for a pure element, with the isotope. By contrast, for X-ray diffraction, the property corresponding to *b* is the atomic scattering factor *A*(*q*), which depends on the direction of observation and electron density in the isolated atom. The structure factor *S*(*q*) obtained by X-ray diffraction has, in general, better accuracy at intermediate values of *q*. At the ends of the scale of *q*, it is less precise than the structure factor obtained by neutron diffraction, because the atomic scattering factor *A*(*q*) is very small for high values of *q* and very poorly known for low values of *q*.

#### **3.2 Structure factor and pair correlation function**

When a photon of wave vector **k** = <sup>2</sup>*<sup>π</sup> <sup>λ</sup>* **u** interacts with an atom, it is deflected by an angle *θ* and the wave vector of the scattered photon is **k**� = <sup>2</sup>*<sup>π</sup> <sup>λ</sup>* **u**� , where **u** and **u**� are unit vectors. If the scattering is elastic it results that |**k**� <sup>|</sup> <sup>=</sup> <sup>|</sup>**k**<sup>|</sup> , because *<sup>E</sup>* <sup>∝</sup> *<sup>k</sup>*<sup>2</sup> <sup>=</sup> *cte*, and that the scattering vector (or transfer vector) **q** is defined by the Bragg law:

$$\mathbf{q} = \mathbf{k'} - \mathbf{k}, \qquad \text{and} \qquad |\mathbf{q}| = 2|\mathbf{k}|\sin\frac{\theta}{2} = \frac{4\pi}{\lambda}\sin\frac{\theta}{2}. \tag{3}$$

6 Thermodynamics book 1

The pair correlation function *g*(*r*) can be deduced from the experimental measurement of the structure factor *S*(*q*) by X-ray, neutron or electron diffractions. In condensed matter, the scatterers are essentially individual atoms, and diffraction experiments can only measure the structure of monatomic liquids such as rare gases and metals. By contrast, they provide no information on the structure of molecular liquids, unless they are composed of spherical

Furthermore, each type of radiation-matter interaction has its own peculiarities. While the electrons are diffracted by all the charges in the atoms (electrons and nuclei), neutrons are diffracted by nuclei and X-rays are diffracted by the electrons localized on stable electron shells. The electron diffraction is practically used for fluids of low density, whereas the beams of neutrons and X-rays are used to study the structure of liquids, with their advantages and disadvantages. For example, the radius of the nuclei being 10, 000 times smaller than that of atoms, it is not surprising that the structure factors obtained with neutrons are not completely

To achieve an experience of X-ray diffraction, we must irradiate the liquid sample with a monochromatic beam of X-rays having a wavelength in the range of the interatomic distance

larger than the mean energy of atoms that is of the order of few *kBT*, namely about 10−<sup>1</sup> eV. The large difference of the masses and energies between a photon and an atom makes that the photon-atom collision is elastic (constant energy) and that the liquid is transparent to the radiation. Naturally, the dimensions of the sample must be sufficiently large compared to the wavelength *λ* of the radiation, so that there are no side effects due to the walls of the enclosure - but not too much though for avoiding excessive absorption of the radiation. This would be particularly troublesome if the X-rays had to pass across metallic elements with large atomic

The incident radiation is characterized by its wavelength *λ* and intensity *I*0, and the diffraction patterns depend on the structural properties of the liquids and on the diffusion properties of atoms. In neutron scattering, the atoms are characterized by the scattering cross section *σ* = <sup>4</sup>*πb*2, where *<sup>b</sup>* is a parameter approximately equal to the radius of the core (<sup>∼</sup> <sup>10</sup>−<sup>14</sup> m). Note that the parameter *b* does not depend on the direction of observation but may vary slightly, even for a pure element, with the isotope. By contrast, for X-ray diffraction, the property corresponding to *b* is the atomic scattering factor *A*(*q*), which depends on the direction of observation and electron density in the isolated atom. The structure factor *S*(*q*) obtained by X-ray diffraction has, in general, better accuracy at intermediate values of *q*. At the ends of the scale of *q*, it is less precise than the structure factor obtained by neutron diffraction, because the atomic scattering factor *A*(*q*) is very small for high values of *q* and very poorly known for

*<sup>λ</sup>* **u** interacts with an atom, it is deflected by an angle *θ*

<sup>|</sup> <sup>=</sup> <sup>|</sup>**k**<sup>|</sup> , because *<sup>E</sup>* <sup>∝</sup> *<sup>k</sup>*<sup>2</sup> <sup>=</sup> *cte*, and that the scattering

*<sup>λ</sup>* sin *<sup>θ</sup>* 2

<sup>2</sup> <sup>=</sup> <sup>4</sup>*<sup>π</sup>*

, where **u** and **u**� are unit vectors. If

. (3)

*<sup>λ</sup>* **u**�

*<sup>λ</sup>* <sup>∼</sup> <sup>10</sup><sup>4</sup> eV), much

(*<sup>λ</sup>* <sup>∼</sup> 0, 1 nm). At this radiation corresponds a photon energy (*h<sup>ν</sup>* <sup>=</sup> *hc*

**3. The structure of liquids**

**3.1 Scattered radiation in liquids**

identical to those obtained with X-rays.

numbers.

low values of *q*.

**3.2 Structure factor and pair correlation function**

and the wave vector of the scattered photon is **k**� = <sup>2</sup>*<sup>π</sup>*

vector (or transfer vector) **q** is defined by the Bragg law:

**<sup>q</sup>** <sup>=</sup> **<sup>k</sup>**� <sup>−</sup> **<sup>k</sup>**, and <sup>|</sup>**q**<sup>|</sup> <sup>=</sup> <sup>2</sup> <sup>|</sup>**k**<sup>|</sup> sin *<sup>θ</sup>*

When a photon of wave vector **k** = <sup>2</sup>*<sup>π</sup>*

the scattering is elastic it results that |**k**�

molecules or monatomic ions, like in some molten salts.

Now, if we consider an assembly of *N* identical atoms forming the liquid sample, the intensity scattered by the atoms in the direction *θ* (or **q**, according to Bragg's law) is given by:

$$I(q) = A\_N A\_N^\star = A\_0 A\_0^\star \sum\_{j=1}^N \sum\_{l=1}^N \exp\left[i\mathbf{q}\left(\mathbf{r}\_j - \mathbf{r}\_l\right)\right].$$

In a crystalline solid, the arrangement of atoms is known once and for all, and the representation of the scattered intensity *I* is given by spots forming the Laue or Debye-Scherrer patterns. But in a liquid, the atoms are in continous motion, and the diffraction experiment gives only the mean value of successive configurations during the experiment. Given the absence of translational symmetry in liquids, this mean value provides no information on long-range order. By contrast, it is a good measure of short-range order around each atom chosen as origin. Thus, in a liquid, the scattered intensity must be expressed as a function of *q* by the statistical average:

$$I(q) = I\_0 \left\langle \sum\_{l=j=1}^{N} \exp\left[i\mathbf{q}\left(\mathbf{r}\_j - \mathbf{r}\_l\right)\right] \right\rangle + I\_0 \left\langle \sum\_{j=1}^{N} \sum\_{l \neq j}^{N} \exp\left[i\mathbf{q}\left(\mathbf{r}\_j - \mathbf{r}\_l\right)\right] \right\rangle. \tag{4}$$

The first mean value, for *l* = *j*, is worth *N* because it represents the sum of *N* terms, each of them being equal to unity. To evaluate the second mean value, one should be able to calculate the sum of exponentials by considering all pairs of atoms (*j*, *l*) in all configurations counted during the experiment, then carry out the average of all configurations. However, this calculation can be achieved only by numerical simulation of a system made of a few particles. In a real system, the method adopted is to determine the mean contribution brought in by each pair of atoms (*j*, *l*), using the probability of finding the atoms *j* and *l* in the positions **r**� and **r**, respectively. To this end, we rewrite the double sum using the Dirac delta function in order to calculate the statistical average in terms of the *density of probability PN*(**r***N*, **p***N*) of the *canonical ensemble*1. Therefore, the statistical average can be written by using the *distribution*

<sup>1</sup> It seems useful to remember that the *probability density* function in the canonical ensemble is:

$$P\_N(\mathbf{r}^N, \mathbf{p}^N) = \frac{1}{N! h^{3N} Q\_N(V, T)} \exp\left[-\beta H\_N(\mathbf{r}^N, \mathbf{p}^N)\right].$$

where *HN*(**r***N*, **<sup>p</sup>***N*) = <sup>∑</sup> *<sup>p</sup>*<sup>2</sup> <sup>2</sup>*<sup>m</sup>* <sup>+</sup> *<sup>U</sup>*(**r***N*) is the *Hamiltonian* of the system, *<sup>β</sup>* <sup>=</sup> <sup>1</sup> *kBT* and *QN*(*V*, *<sup>T</sup>*) the *partition function* defined as:

$$Q\_N(V,T) = \frac{Z\_N(V,T)}{N!\Lambda^{3N}}\prime$$

with the *thermal wavelength* Λ, which is a measure of the thermodynamic uncertainty in the localization of a particle of mass *m*, and the *configuration integral ZN*(*V*, *T*), which is expressed in terms of the total *potential energy U*(**r***N*). They read:

$$\begin{aligned} \Lambda &= \sqrt{\frac{\hbar^2}{2\pi m k\_B T}}, \\ \text{and} \quad Z\_N(V, T) &= \int \int\_N \exp\left[-\beta \mathcal{U}(\mathbf{r}^N)\right] d\mathbf{r}^N. \end{aligned}$$

Besides, the partition function *QN*(*V*, *T*) allows us to determine the free energy *F* according to the relation:

$$F = E - TS = -k\_B T \ln Q\_N(V, T).$$

The reader is advised to consult statistical-physics textbooks for further details.

of Simple Liquids 9

Thermodynamic Perturbation Theory of Simple Liquids 847

Fig. 3. Structure factor *S*(*q*) and pair correlation function *g*(*r*) of simple liquids.

function, the expression of the scattered intensity *I*(*q*) becomes:

 *V*

*structure factor S*(*q*) is defined by the following normalized function:

*I*(*q*) = *N I*<sup>0</sup> + *N I*0*ρ*

*<sup>S</sup>*(*q*) = *<sup>I</sup>*(*q*) <sup>−</sup> (2*π*)3*N I*0*ρδ*(**q**) *N I*<sup>0</sup>

exact relation *S*(0) = *ρkBTχT*.

which is zero for all values of *q*, except in *q* = 0 for which it is infinite. In using the delta

From the experimental point of view, it is necessary to exclude the measurement of the scattered intensity in the direction of the incident beam (*q* = 0). Therefore, in practice, the

Consequently, the pair correlation function *g*(*r*) can be extracted from the experimental results

*ρ* [*g*(*r*) − 1] = *TF* [*S*(*q*) − 1] . The pair correlation function *g*(*r*) is a dimensionless quantity, whose the graphic representation is given in figure (3). The gap around unity measures the probability of finding a particle at distance *r* from a particle taken in an arbitrary origin. The main peak of *g*(*r*) corresponds to the position of first neighbors, and the successive peaks to the next close neighbors. The pair correlation function *g*(*r*) clearly shows the existence of a short-range order that is fading rapidly beyond four or five interatomic distances. In passing, it should be mentioned that the structure factor at *q* = 0 is related to the isothermal compressibility by the

of the structure factor *S*(*q*) by performing the numerical Fourier transformation:

<sup>=</sup> <sup>1</sup> <sup>+</sup> <sup>4</sup>*πρ* <sup>∞</sup>

0

sin(*qr*)

*qr* [*g*(*r*) <sup>−</sup> <sup>1</sup>]*<sup>r</sup>*

<sup>2</sup>*dr*. (7)

exp (*i***qr**) [*g*(*r*) <sup>−</sup> <sup>1</sup>] *<sup>d</sup>***<sup>r</sup>** <sup>+</sup> *N I*0*ρ*(2*π*)3*δ*(**q**).

*function*<sup>2</sup> *ρ* (2) *<sup>N</sup>* (**r**,**r**� ) in the form:

$$\left\langle \sum\_{j=1}^{N} \sum\_{l \neq j}^{N} \exp\left[i\mathbf{q}\left(\mathbf{r}\_{j} - \mathbf{r}\_{l}\right)\right] \right\rangle = \int \int\_{6} d\mathbf{r} d\mathbf{r}' \exp\left[i\mathbf{q}\left(\mathbf{r}' - \mathbf{r}\right)\right] \rho\_{N}^{(2)}\left(\mathbf{r}, \mathbf{r}'\right) \dots$$

If the liquid is assumed to be homogeneous and isotropic, and that all atoms have the same properties, one can make the changes of variables **R** = **r** and **X** = **r**� − **r**, and explicit the pair correlation function *<sup>g</sup>*(|**r**� <sup>−</sup> **<sup>r</sup>**|) = *<sup>ρ</sup>* (2) *<sup>N</sup>* (**r**,**r**� ) *<sup>ρ</sup>*<sup>2</sup> in the statistical average as3:

$$\left\langle \sum\_{j=1}^{N} \sum\_{l \neq j}^{N} \exp\left[i\mathbf{q}\left(\mathbf{r}\_{j} - \mathbf{r}\_{l}\right)\right] \right\rangle = 4\pi\rho^{2}V \int\_{0}^{\infty} \frac{\sin(qr)}{qr} \mathbf{g}(r)r^{2} dr. \tag{5}$$

One sees that the previous integral diverges because the integrand increases with *r*. The problem comes from the fact that the scattered intensity, for *q* = 0, has no physical meaning and can not be measured. To overcome this difficulty, one rewrites the scattered intensity *I*(*q*) defined by equation (4) in the equivalent form (cf. footnote 3):

$$I(q) = Nl\_0 + Nl\_0 \rho \int\_V \exp\left(i\mathbf{q}\mathbf{r}\right) \left[\mathbf{g}(r) - 1\right] d\mathbf{r} + Nl\_0 \rho \int\_V \exp\left(i\mathbf{q}\mathbf{r}\right) d\mathbf{r}.\tag{6}$$

To large distances, *g*(*r*) tends to unity, so that [*g*(*r*) − 1] tends towards zero, making the first integral convergent. As for the second integral, it corresponds to the Dirac delta function4,

<sup>2</sup> It should be stressed that the distribution function *ρ* (2) *N* **r**2 is expressed as:

$$\rho\_N^{(2)}\left(\mathbf{r},\mathbf{r}'\right) = \rho^2 \mathcal{g}(|\mathbf{r}'-\mathbf{r}|) = \frac{N!}{(N-2)!Z\_N} \int \int\_{\mathfrak{J}(N-2)} \exp\left[-\beta \mathcal{U}(\mathbf{r}^N)\right] d\mathbf{r}\_3...d\mathbf{r}\_N.$$

<sup>3</sup> To evaluate an integral of the form:

$$I = \int\_{V} d\mathbf{r} \exp\left(i\mathbf{q}\mathbf{r}\right) \mathcal{J}(r),$$

one must use the spherical coordinates by placing the vector **q** along the *z* axis, where *θ* = (**q**,**r**). Thus, the integral reads:

$$I = \int\_0^{2\pi} d\varphi \int\_0^{\pi} \int\_0^{\infty} \exp\left(iqr\cos\theta\right) \mathcal{g}(r)r^2 \sin\theta d\theta dr$$

with *μ* = cos *θ* and *dμ* = − sin *θdθ*. It follows that:

$$I = -2\pi \int\_0^\infty \left[ \int\_{+1}^{-1} \exp\left(iqr\mu\right) d\mu \right] \mathcal{g}(r) r^2 dr = 4\pi \int\_0^\infty \frac{\sin(qr)}{qr} \mathcal{g}(r) r^2 dr.$$

<sup>4</sup> The generalization of the Fourier transform of the Dirac delta function to three dimensions is:

$$\delta(\mathbf{r}) = \frac{1}{\left(2\pi\right)^{3}} \int \int \int\_{-\infty}^{+\infty} \delta(\mathbf{q}) \exp\left(-i\mathbf{q}\mathbf{r}\right) d\mathbf{q} = \frac{1}{\left(2\pi\right)^{3}}.$$

and the inverse transform is:

$$\delta(\mathbf{q}) = \int \int \int\_{-\infty}^{+\infty} \delta(\mathbf{r}) \exp\left(i\mathbf{q}\mathbf{r}\right) d\mathbf{r} = \frac{1}{\left(2\pi\right)^{3}} \int \int \int\_{-\infty}^{+\infty} \exp\left(i\mathbf{q}\mathbf{r}\right) d\mathbf{r} \dots$$

8 Thermodynamics book 1

If the liquid is assumed to be homogeneous and isotropic, and that all atoms have the same properties, one can make the changes of variables **R** = **r** and **X** = **r**� − **r**, and explicit the pair

One sees that the previous integral diverges because the integrand increases with *r*. The problem comes from the fact that the scattered intensity, for *q* = 0, has no physical meaning and can not be measured. To overcome this difficulty, one rewrites the scattered intensity *I*(*q*)

To large distances, *g*(*r*) tends to unity, so that [*g*(*r*) − 1] tends towards zero, making the first integral convergent. As for the second integral, it corresponds to the Dirac delta function4,

> ) = *<sup>N</sup>*! (*N* − 2)!*ZN*

> > *I* = *V*

> > > ∞ 0

exp (*iqrμ*) *dμ*

 +∞ −∞

<sup>4</sup> The generalization of the Fourier transform of the Dirac delta function to three dimensions is:

*<sup>δ</sup>*(**r**) exp (*i***qr**) *<sup>d</sup>***<sup>r</sup>** <sup>=</sup> <sup>1</sup>

� exp *i***q r** � − **r** *ρ* (2) *N* **r**,**r** � .

*<sup>ρ</sup>*<sup>2</sup> in the statistical average as3:

 ∞ 0

sin(*qr*) *qr <sup>g</sup>*(*r*)*<sup>r</sup>*

> *V*

is expressed as:

exp 

<sup>−</sup>*βU*(**r***N*)

<sup>2</sup> sin *θdθdr*,

sin(*qr*) *qr <sup>g</sup>*(*r*)*<sup>r</sup>*

> (2*π*) 3 ,

> > exp (*i***qr**) *d***r**.

 ∞ 0

 +∞ −∞

*d***r**3...*d***r***N*.

<sup>2</sup>*dr*.

<sup>2</sup>*dr*. (5)

exp (*i***qr**) *d***r**. (6)

= 4*πρ*2*V*

exp (*i***qr**) [*g*(*r*) − 1] *d***r** + *N I*0*ρ*

(2) *N* **r**2

*d***r** exp (*i***qr**) *g*(*r*),

exp (*iqr* cos *θ*) *g*(*r*)*r*

<sup>2</sup>*dr* = 4*π*

*<sup>δ</sup>*(**q**) exp (−*i***qr**) *<sup>d</sup>***<sup>q</sup>** <sup>=</sup> <sup>1</sup>

(2*π*) 3

one must use the spherical coordinates by placing the vector **q** along the *z* axis, where *θ* = (**q**,**r**). Thus,

 *g*(*r*)*r* 3(*N*−2)

 = 6 *d***r***d***r**

(2) *<sup>N</sup>* (**r**,**r**� )

defined by equation (4) in the equivalent form (cf. footnote 3):

 *V*

*function*<sup>2</sup> *ρ*

(2) *<sup>N</sup>* (**r**,**r**�

> *N* ∑ *j*=1

*N* ∑ *l*�=*j* exp *i***q r***<sup>j</sup>* − **r***<sup>l</sup>*

correlation function *<sup>g</sup>*(|**r**� <sup>−</sup> **<sup>r</sup>**|) = *<sup>ρ</sup>*

 *N* ∑ *j*=1

*N* ∑ *l*�=*j* exp *i***q r***<sup>j</sup>* − **r***<sup>l</sup>*

*I*(*q*) = *N I*<sup>0</sup> + *N I*0*ρ*

<sup>2</sup> It should be stressed that the distribution function *ρ*

*I* = 2*π* 0 *dϕ π* 0

 ∞ 0

*<sup>δ</sup>*(**r**) = <sup>1</sup> (2*π*) 3

 +∞ −∞

 −<sup>1</sup> +1

with *μ* = cos *θ* and *dμ* = − sin *θdθ*. It follows that:

*I* = −2*π*

*δ*(**q**) =

and the inverse transform is:

*ρ* (2) *N* **r**,**r** � = *ρ*<sup>2</sup> *g*( **r** � − **r** 

the integral reads:

<sup>3</sup> To evaluate an integral of the form:

) in the form:

Fig. 3. Structure factor *S*(*q*) and pair correlation function *g*(*r*) of simple liquids.

which is zero for all values of *q*, except in *q* = 0 for which it is infinite. In using the delta function, the expression of the scattered intensity *I*(*q*) becomes:

$$I(q) = NI\_0 + NI\_0 \rho \int\_V \exp\left(i\mathbf{q}\mathbf{r}\right) \left[\mathbf{g}(r) - \mathbf{1}\right] d\mathbf{r} + NI\_0 \rho (2\pi)^3 \delta(\mathbf{q}).$$

From the experimental point of view, it is necessary to exclude the measurement of the scattered intensity in the direction of the incident beam (*q* = 0). Therefore, in practice, the *structure factor S*(*q*) is defined by the following normalized function:

$$S(q) = \frac{I(q) - (2\pi)^3 N I\_0 \rho \delta(\mathbf{q})}{N I\_0} = 1 + 4\pi \rho \int\_0^\infty \frac{\sin(qr)}{qr} \left[g(r) - 1\right] r^2 dr. \tag{7}$$

Consequently, the pair correlation function *g*(*r*) can be extracted from the experimental results of the structure factor *S*(*q*) by performing the numerical Fourier transformation:

$$
\rho \left[ \mathbf{g}(r) - 1 \right] = T F \left[ \mathbf{S}(q) - 1 \right].
$$

The pair correlation function *g*(*r*) is a dimensionless quantity, whose the graphic representation is given in figure (3). The gap around unity measures the probability of finding a particle at distance *r* from a particle taken in an arbitrary origin. The main peak of *g*(*r*) corresponds to the position of first neighbors, and the successive peaks to the next close neighbors. The pair correlation function *g*(*r*) clearly shows the existence of a short-range order that is fading rapidly beyond four or five interatomic distances. In passing, it should be mentioned that the structure factor at *q* = 0 is related to the isothermal compressibility by the exact relation *S*(0) = *ρkBTχT*.

of Simple Liquids 11

Thermodynamic Perturbation Theory of Simple Liquids 849

Therefore, the calculation of internal energy of a liquid requires knowledge of the pair potential *u*(*r*) and the pair correlation function *g*(*r*). For the latter, the choice is to employ either the experimental values or values derived from the microscopic theory of liquids. Note that the integrand in equation (9) is the product of the pair potential by the pair correlation function, weighted by *r*2. It should be also noted that the calculation of *E* can be made taking into account the three-body potential *u*3(**r**1,**r**2,**r**3) and the three-body correlation function *g*(3)(**r**1,**r**2,**r**3). In this case, the correlation function at three bodies must be determined only

The expression of the pressure is obtained in the same way that the internal energy, in

The derivation of the configuration integral with respect to volume requires using the reduced

volume. Indeed, if the volume element is *d***r** = *Vd***X**, the scalar variable *dr* = *V*1/3*dX* leads to

In view of this, the configuration integral and its derivative with respect to *V* are written in

Assuming that the potential energy is decomposed into a sum of pair potentials, and with the help of equation (10), the derivation of the potential energy versus volume is performed as:

> <sup>3</sup>*<sup>V</sup>* ∑ *i* ∑ *j*>*i rij*

> > *rij*

Like for the calculation of internal energy, the additivity assumption of pair potentials permits

<sup>−</sup>*βU*(**r***N*)

 3*N*

exp

<sup>−</sup>*βU*(**r***N*)

 3*N* 

*<sup>V</sup>*−2/3*<sup>X</sup>* <sup>=</sup> <sup>1</sup>

*d***X**1...*d***X***N*,

 exp

*∂u*(*rij*) *∂rij* ,

exp

*<sup>d</sup>***r***<sup>N</sup>* <sup>=</sup> *<sup>N</sup>*(*<sup>N</sup>* <sup>−</sup> <sup>1</sup>) 2

*∂u*(*rij*) *∂rij*

<sup>−</sup>*<sup>β</sup> <sup>∂</sup>U*(**r***N*) *∂V*

*<sup>V</sup>* that allows us to find the dependence of the potential energy *<sup>U</sup>*(**r***N*) versus

*<sup>∂</sup><sup>V</sup>* ln *ZN*(*V*, *<sup>T</sup>*).

<sup>3</sup>*<sup>V</sup> <sup>r</sup>*. (10)

<sup>−</sup>*βU*(**r***N*)

<sup>−</sup>*βU*(**r***N*)

 *rij*

*d***X**1...*d***X***N*.

*d***r**1...*d***r***N*. (11)

*∂u*(*rij*) *<sup>∂</sup>rij* ,

*<sup>∂</sup><sup>V</sup>* ln *QN*(*V*, *<sup>T</sup>*) = *kBT <sup>∂</sup>*

by the theory of liquids (7), since it is not accessible by experiment.

*dr dV* <sup>=</sup> <sup>1</sup> 3

exp

*V<sup>N</sup> ZN*(*V*, *T*)

> *∂U*(**r***N*) *<sup>∂</sup><sup>V</sup>* <sup>=</sup> <sup>1</sup>

> > *i* ∑ *j*>*i*

*∂u*(*rij*) *∂rij*

*<sup>p</sup>* <sup>=</sup> *kBT <sup>∂</sup>*

 3*N*

*<sup>V</sup>* <sup>+</sup>

so that the expression of the pressure becomes:

1 *ZN*(*V*, *<sup>T</sup>*) ∑

us to write the sum of integrals of the previous equation as:

*rij*

*<sup>V</sup>* <sup>−</sup> <sup>1</sup> 3*V*

1 *ZN*(*V*, *T*)

the following forms with reduced variables:

*ZN*(*V*, *T*) = *V<sup>N</sup>*

*<sup>∂</sup><sup>V</sup>* ln *ZN*(*V*, *<sup>T</sup>*) = *<sup>N</sup>*

*<sup>p</sup>* <sup>=</sup> *kBT <sup>N</sup>*

∑ *i* ∑ *j*>*i*

**4.2 Pressure**

variable **X** = **<sup>r</sup>**

the derivative:

*∂*

considering the equation:

### **4. Thermodynamic functions of liquids**

#### **4.1 Internal energy**

To express the internal energy of a liquid in terms of the pair correlation function, one must first use the following relation from statistical mechanics :

$$E = k\_B T^2 \frac{\partial}{\partial T} \ln Q\_N(V, T)\_V$$

where the partition function *QN*(*V*, *T*) depends on the configuration integral *ZN*(*V*, *T*) and on the thermal wavelength Λ, in accordance with the equations given in footnote (1). The derivative of ln *QN*(*V*, *T*) with respect to *T* can be written:

$$\frac{\partial}{\partial T} \ln Q\_N(V, T) = \frac{\partial}{\partial T} \ln Z\_N(V, T) - 3N \frac{\partial}{\partial T} \ln \Lambda\_\prime$$

with:

$$\begin{split} \frac{\partial}{\partial T} \ln Z\_N(V, T) &= \frac{1}{Z\_N(V, T)} \int \int \left[ \frac{1}{k\_B T^2} U(\mathbf{r}^N) \right] \exp\left[ -\beta U(\mathbf{r}^N) \right] d\mathbf{r}^N \\ \text{and} \qquad \frac{\partial}{\partial T} \ln \Lambda &= \frac{1}{\Lambda} \left( -\frac{1}{2T^{3/2}} \sqrt{\frac{h^2}{2\pi m k\_B}} \right) = -\frac{1}{2T} .\end{split}$$

Then, the calculation is continued by admitting that the total potential energy *U*(**r***N*) is written as a sum of pair potentials, in the form *U*(**r***N*) = ∑*<sup>i</sup>* ∑*j*>*<sup>i</sup> u*(*rij*). The internal energy reads:

$$E = \frac{3}{2} \mathcal{N} k\_B T + \frac{1}{\mathcal{Z}\_N(V, T)} \int \int \left[ \sum\_i \sum\_{j>i} u(r\_{ij}) \right] \exp\left[ -\beta \mathcal{U}(\mathbf{r}^N) \right] d\mathbf{r}^N. \tag{8}$$

The first term on the RHS corresponds to the kinetic energy of the system; it is the ideal gas contribution. The second term represents the potential energy. Given the assumption of additivity of pair potentials, we can assume that it is composed of *N*(*N* − 1)/2 identical terms, permitting us to write:

$$\sum\_{i} \sum\_{j>i} \frac{1}{Z\_N(V, T)} \int \int u(r\_{ij}) \exp\left[-\beta U(\mathbf{r}^N)\right] d\mathbf{r}^N = \frac{N(N-1)}{2} \left< u(r\_{ij}) \right> ,$$

where the mean value is expressed in terms of the pair correlation function as:

$$
\langle \mu(r\_{12}) \rangle = \rho^2 \frac{(N-2)!}{N!} \int \int\_{\mathfrak{G}} \mu(r\_{12}) \left[ \mathfrak{g}\_N^{(2)}(\mathbf{r}\_1, \mathbf{r}\_2) \right] d\mathbf{r}\_1 d\mathbf{r}\_2.
$$

For a homogeneous and isotropic fluid, one can perform the change of variables **R** = **r**<sup>1</sup> and **r** = **r**<sup>1</sup> − **r**2, where **R** and **r** describe the system volume, and write the expression of internal energy in the integral form:

$$E = \frac{3}{2} N k\_B T + 2\pi \rho N \int\_0^\infty u(r) g(r) r^2 dr. \tag{9}$$

Therefore, the calculation of internal energy of a liquid requires knowledge of the pair potential *u*(*r*) and the pair correlation function *g*(*r*). For the latter, the choice is to employ either the experimental values or values derived from the microscopic theory of liquids. Note that the integrand in equation (9) is the product of the pair potential by the pair correlation function, weighted by *r*2. It should be also noted that the calculation of *E* can be made taking into account the three-body potential *u*3(**r**1,**r**2,**r**3) and the three-body correlation function *g*(3)(**r**1,**r**2,**r**3). In this case, the correlation function at three bodies must be determined only by the theory of liquids (7), since it is not accessible by experiment.

#### **4.2 Pressure**

10 Thermodynamics book 1

To express the internal energy of a liquid in terms of the pair correlation function, one must

where the partition function *QN*(*V*, *T*) depends on the configuration integral *ZN*(*V*, *T*) and on the thermal wavelength Λ, in accordance with the equations given in footnote (1). The

*<sup>∂</sup><sup>T</sup>* ln *QN*(*V*, *<sup>T</sup>*),

*<sup>∂</sup><sup>T</sup>* ln *ZN*(*V*, *<sup>T</sup>*) <sup>−</sup> <sup>3</sup>*<sup>N</sup> <sup>∂</sup>*

*h*2 2*πmkB*

*u*(*rij*)

�

*<sup>d</sup>***r***<sup>N</sup>* <sup>=</sup> *<sup>N</sup>*(*<sup>N</sup>* <sup>−</sup> <sup>1</sup>) 2

> � *d***r**1*d***r**2.

*kBT*<sup>2</sup> *<sup>U</sup>*(**r***N*)

⎞ <sup>⎠</sup> <sup>=</sup> <sup>−</sup> <sup>1</sup> 2*T* .

⎤ ⎦ exp �

� � � 1

�

Then, the calculation is continued by admitting that the total potential energy *U*(**r***N*) is written as a sum of pair potentials, in the form *U*(**r***N*) = ∑*<sup>i</sup>* ∑*j*>*<sup>i</sup> u*(*rij*). The internal energy reads:

The first term on the RHS corresponds to the kinetic energy of the system; it is the ideal gas contribution. The second term represents the potential energy. Given the assumption of additivity of pair potentials, we can assume that it is composed of *N*(*N* − 1)/2 identical terms,

<sup>−</sup>*βU*(**r***N*)

� � <sup>⎡</sup> <sup>⎣</sup>∑ *i* ∑ *j*>*i*

�

� � 6 *u*(*r*12) � *g* (2 *<sup>N</sup>* (**r**1,**r**2)

For a homogeneous and isotropic fluid, one can perform the change of variables **R** = **r**<sup>1</sup> and **r** = **r**<sup>1</sup> − **r**2, where **R** and **r** describe the system volume, and write the expression of internal

> � ∞ 0

*u*(*r*)*g*(*r*)*r*

*<sup>∂</sup><sup>T</sup>* ln <sup>Λ</sup>,

<sup>−</sup>*βU*(**r***N*)

<sup>−</sup>*βU*(**r***N*)

�

� *u*(*rij*) � ,

<sup>2</sup>*dr*. (9)

*d***r***N*. (8)

� *d***r***<sup>N</sup>*

� exp �

*<sup>E</sup>* <sup>=</sup> *kBT*<sup>2</sup> *<sup>∂</sup>*

*ZN*(*V*, *T*)

⎛ <sup>⎝</sup><sup>−</sup> <sup>1</sup> 2*T*3/2

Λ

1 *ZN*(*V*, *T*)

**4. Thermodynamic functions of liquids**

first use the following relation from statistical mechanics :

derivative of ln *QN*(*V*, *T*) with respect to *T* can be written:

*<sup>∂</sup><sup>T</sup>* ln *ZN*(*V*, *<sup>T</sup>*) = <sup>1</sup>

*<sup>∂</sup><sup>T</sup>* ln <sup>Λ</sup> <sup>=</sup> <sup>1</sup>

*<sup>∂</sup><sup>T</sup>* ln *QN*(*V*, *<sup>T</sup>*) = *<sup>∂</sup>*

*∂*

*∂*

and *<sup>∂</sup>*

*<sup>E</sup>* <sup>=</sup> <sup>3</sup> 2

permitting us to write:

∑ *i* ∑ *j*>*i*

energy in the integral form:

*NkBT* +

1 *ZN*(*V*, *T*) � �

�*u*(*r*12)� <sup>=</sup> *<sup>ρ</sup>*<sup>2</sup> (*<sup>N</sup>* <sup>−</sup> <sup>2</sup>)!

*<sup>E</sup>* <sup>=</sup> <sup>3</sup> 2

*u*(*rij*) exp

where the mean value is expressed in terms of the pair correlation function as:

*N*!

*NkBT* + 2*πρN*

**4.1 Internal energy**

with:

The expression of the pressure is obtained in the same way that the internal energy, in considering the equation:

$$p = k\_B T \frac{\partial}{\partial V} \ln Q\_N(V, T) = k\_B T \frac{\partial}{\partial V} \ln Z\_N(V, T).$$

The derivation of the configuration integral with respect to volume requires using the reduced variable **X** = **<sup>r</sup>** *<sup>V</sup>* that allows us to find the dependence of the potential energy *<sup>U</sup>*(**r***N*) versus volume. Indeed, if the volume element is *d***r** = *Vd***X**, the scalar variable *dr* = *V*1/3*dX* leads to the derivative:

$$\frac{dr}{dV} = \frac{1}{3}V^{-2/3}X = \frac{1}{3V}r.\tag{10}$$

In view of this, the configuration integral and its derivative with respect to *V* are written in the following forms with reduced variables:

$$\begin{split} Z\_{N}(V,T) &= V^{N} \int \int\_{3N} \exp\left[-\beta \mathcal{U}(\mathbf{r}^{N})\right] d\mathbf{X}\_{1}...d\mathbf{X}\_{N}, \\ \frac{\partial}{\partial V} \ln Z\_{N}(V,T) &= \frac{N}{V} + \frac{V^{N}}{Z\_{N}(V,T)} \int \int\_{3N} \left[-\beta \frac{\partial \mathcal{U}(\mathbf{r}^{N})}{\partial V}\right] \exp\left[-\beta \mathcal{U}(\mathbf{r}^{N})\right] d\mathbf{X}\_{1}...d\mathbf{X}\_{N}. \end{split}$$

Assuming that the potential energy is decomposed into a sum of pair potentials, and with the help of equation (10), the derivation of the potential energy versus volume is performed as:

$$\frac{\partial \mathcal{U}(\mathbf{r}^N)}{\partial V} = \frac{1}{\Im V} \sum\_{\mathbf{i}} \sum\_{j>\mathbf{i}} r\_{\mathbf{ij}} \frac{\partial u(r\_{\mathbf{ij}})}{\partial r\_{\mathbf{ij}}} \lambda$$

so that the expression of the pressure becomes:

$$p = k\_B T \frac{N}{V} - \frac{1}{3V} \frac{1}{Z\_N(V, T)} \sum\_i \sum\_{j>i} \int \int\_{3N} \left[ r\_{ij} \frac{\partial u(r\_{ij})}{\partial r\_{ij}} \right] \exp\left[ -\beta l I(\mathbf{r}^N) \right] d\mathbf{r}\_1...d\mathbf{r}\_N. \tag{11}$$

Like for the calculation of internal energy, the additivity assumption of pair potentials permits us to write the sum of integrals of the previous equation as:

$$\sum\_{i} \sum\_{j>i} \frac{1}{Z\_N(V, T)} \int \int \left[ r\_{ij} \frac{\partial u(r\_{ij})}{\partial r\_{ij}} \right] \exp\left[ -\beta \mathcal{U}(\mathbf{r}^N) \right] d\mathbf{r}^N = \frac{N(N-1)}{2} \left\langle r\_{ij} \frac{\partial u(r\_{ij})}{\partial r\_{ij}} \right\rangle\_{\mathbf{r}^N}$$

of Simple Liquids 13

Thermodynamic Perturbation Theory of Simple Liquids 851

means that particle 1 is completely coupled with the other particles, while *λ* = 0 indicates a zero coupling, that is to say the absence of the particle 1 in the system. This allows the writing

*u*(*rij*) =

*N* ∑ *i*

*u*(*rij*) = *U*(**r***N*−1).

<sup>−</sup>*βU*(**r***N*−1)

*N* ∑ *j*>*i*≥1

*u*(*rij*) = *U*(**r***N*),

*d***r**1*d***r**2...*d***r***<sup>N</sup>* = *ZN*(*V*, *T*), (15)

*ZN*(*V*, *<sup>T</sup>*, *<sup>λ</sup>* <sup>=</sup> <sup>0</sup>) <sup>+</sup> ln *<sup>V</sup>* (17)

*<sup>N</sup>* (**r**1,**r**2, *λ*)

<sup>2</sup>*dr*,

*∂λ <sup>d</sup>λ*. (18)

*∂λ* . In particular, with the result

<sup>2</sup>*drdλ*. (19)

*d***r**1*d***r**2.

*<sup>d</sup>***r**2...*d***r***<sup>N</sup>* = *VZN*−1(*V*, *<sup>T</sup>*). (16)

of the important relations:

and

*U*(**r***N*, 1) =

coupling (*λ* = 0) are respectively:

*ZN*(*V*, *<sup>T</sup>*, *<sup>λ</sup>* <sup>=</sup> <sup>1</sup>) =

*ZN*(*V*, *<sup>T</sup>*, *<sup>λ</sup>* <sup>=</sup> <sup>0</sup>) =

integrals in equation (13):

*N* ∑ *j*=2

*u*(*r*1*j*) +

*U*(**r***N*, 0) =

exp

3*N*

ln *ZN*(*V*, *<sup>T</sup>*)

*∂λ* <sup>=</sup> <sup>−</sup>*βρ*<sup>2</sup> (*<sup>N</sup>* <sup>−</sup> <sup>1</sup>)(*<sup>N</sup>* <sup>−</sup> <sup>2</sup>)!

*∂* ln *ZN*(*V*, *T*, *λ*)

arrives to the expression of the chemical potential:

*V d***r**<sup>1</sup> 

(14), we can easily evaluate the partial derivatives *<sup>∂</sup>ZN*

of the footnote (1), we can write *<sup>∂</sup>* ln *ZN*

*∂* ln *ZN*(*V*, *T*, *λ*)

following form:

*N* ∑ *i*

> *N* ∑ *i*

<sup>−</sup>*βU*(**r***N*)

3(*N*−1)

*N* ∑ *j*>*i*≥2

> *N* ∑ *j*>*i*≥2

Under these conditions, the configuration integrals for a total coupling (*λ* = 1) and a zero

exp

*ZN*−1(*V*, *<sup>T</sup>*) <sup>=</sup> ln *ZN*(*V*, *<sup>T</sup>*, *<sup>λ</sup>* <sup>=</sup> <sup>1</sup>)

= ln *V* +

*N*!

*∂λ* <sup>=</sup> <sup>−</sup> *βρ*<sup>2</sup>

*<sup>μ</sup>* <sup>=</sup> *kBT* ln *<sup>ρ</sup>*Λ<sup>3</sup> <sup>+</sup> <sup>4</sup>*πρ* <sup>1</sup>

But with the configuration integral *ZN*(*V*, *T*, *λ*), in which potential energy is given by equation

In addition, if the fluid is homogeneous and isotropic, the above relation simplifies under the

*<sup>N</sup> <sup>V</sup>*

that remains only to be substituted in equation (18) for obtaining the logarithm of the ratio of configuration integrals. And by putting the last expression in equation (13), one ultimately

0

 ∞ 0

 1 0

> 6 *u*(*r*12) *g* (2)

 ∞ 0

*∂* ln *ZN*

*∂λ* and *<sup>∂</sup>* ln *ZN*

*∂λ* as a function of the pair correlation function as:

*u*(*r*)*g*(*r*, *λ*)4*πr*

*u*(*r*)*g*(*r*, *λ*)*r*

These expressions are then used to calculate the logarithm of the ratio of configuration

where the mean value is expressed with the pair correlation function by:

$$\left\langle r\_{12} \frac{\partial u(r\_{12})}{\partial r\_{12}} \right\rangle = \rho^2 \frac{(N-2)!}{N!} \int \int\_{\mathcal{S}} r\_{12} \frac{\partial u(r\_{12})}{\partial r\_{12}} \left[ \mathcal{g}\_N^{(2)}(\mathbf{r}\_1, \mathbf{r}\_2) \right] d\mathbf{r}\_1 d\mathbf{r}\_2.$$

For a homogeneous and isotropic fluid, one can perform the change of variables **R** = **r**<sup>1</sup> and **r** = **r**<sup>1</sup> − **r**2, and simplify the expression of pressure as:

$$p = k\_B T \frac{N}{V} - \frac{2\pi}{3} \rho^2 \int\_0^\infty r^3 \frac{\partial u(r)}{\partial r} g(r) dr. \tag{12}$$

The previous equation provides the pressure of a liquid as a function of the pair potential and the pair correlation function. It is the so-called *pressure equation of state* of liquids. It should be stressed that this equation of state is not unique, as we will see in presenting the hard-sphere reference system (§ 4. 4). As the internal energy, the pressure can be written with an additional term containing the three-body potential *u*3(**r**1,**r**2,**r**3) and the three-body correlation function *g*(3)(**r**1,**r**2,**r**3).

#### **4.3 Chemical potential and entropy**

We are now able to calculate the internal energy (Eq. 9) and pressure (Eq. 12) for any system, of which the potential energy is made of a sum of pair potentials *u*(*r*) and the pair correlation function *g*(*r*) is known. Beside this, all other thermodynamic properties can be easily derived. Traditionally, it is appropriate to derive the chemical potential *μ* as a function of *g*(*r*) by integrating the partition function with respect to a parameter *λ* to be defined (8).

Firstly, the formal expression of the chemical potential is defined by the energy required to introduce a new particle in the system:

$$\mu = F(V, T, N) - F(V, T, N - 1) = \left(\frac{\partial F}{\partial N}\right)\_{V, T}.$$

From footnote (1), the free energy *F* is written:

$$F(V, T, N) = -k\_B T \ln Q\_N(V, T) = -k\_B T \left[ \ln Z\_N(V, T) - \ln N! - N \ln \Lambda^3 \right]\_{\text{max}}$$

so that the chemical potential can be simplified as:

$$\mu = k\_B T \left[ -\ln \frac{Z\_N(V, T)}{Z\_{N-1}(V, T)} + \ln N + \ln \Lambda^3 \right].\tag{13}$$

Secondly, the procedure requires to write the potential energy as a function of the *coupling parameter λ*, under the following form, in order to assess the argument of the logarithm in the above relation:

$$\mathcal{U}(\mathbf{r}^N, \boldsymbol{\lambda}) = \lambda \sum\_{j=2}^N \boldsymbol{\mu}(r\_{1j}) + \sum\_{i}^N \sum\_{j>i \ge 2}^N \boldsymbol{\mu}(r\_{ij}). \tag{14}$$

Varying from 0 to 1, the coupling parameter *λ* measures the degree of coupling of the particle to which it is assigned (1 in this case) with the rest of the system. In the previous relation, *λ* = 1 means that particle 1 is completely coupled with the other particles, while *λ* = 0 indicates a zero coupling, that is to say the absence of the particle 1 in the system. This allows the writing of the important relations:

$$\mathcal{U}(\mathbf{r}^N, 1) = \sum\_{j=2}^N \mu(r\_{1j}) + \sum\_{i}^N \sum\_{j>i \ge 2}^N \mu(r\_{i\bar{j}}) = \sum\_{i}^N \sum\_{j>i \ge 1}^N \mu(r\_{i\bar{j}}) = \mathcal{U}(\mathbf{r}^N).$$

and

12 Thermodynamics book 1

 6 *r*<sup>12</sup>

For a homogeneous and isotropic fluid, one can perform the change of variables **R** = **r**<sup>1</sup> and

The previous equation provides the pressure of a liquid as a function of the pair potential and the pair correlation function. It is the so-called *pressure equation of state* of liquids. It should be stressed that this equation of state is not unique, as we will see in presenting the hard-sphere reference system (§ 4. 4). As the internal energy, the pressure can be written with an additional term containing the three-body potential *u*3(**r**1,**r**2,**r**3) and the three-body correlation function

We are now able to calculate the internal energy (Eq. 9) and pressure (Eq. 12) for any system, of which the potential energy is made of a sum of pair potentials *u*(*r*) and the pair correlation function *g*(*r*) is known. Beside this, all other thermodynamic properties can be easily derived. Traditionally, it is appropriate to derive the chemical potential *μ* as a function of *g*(*r*) by

Firstly, the formal expression of the chemical potential is defined by the energy required to

*ZN*−1(*V*, *<sup>T</sup>*) <sup>+</sup> ln *<sup>N</sup>* <sup>+</sup> ln <sup>Λ</sup><sup>3</sup>

*N* ∑ *i*

*N* ∑ *j*>*i*≥2

 ∞ 0 *r* <sup>3</sup> *∂u*(*r*)

*∂u*(*r*12) *∂r*<sup>12</sup>

 *g* (2 *<sup>N</sup>* (**r**1,**r**2)

 *∂F ∂N* 

*V*,*T* .

ln *ZN*(*V*, *<sup>T</sup>*) <sup>−</sup> ln *<sup>N</sup>*! <sup>−</sup> *<sup>N</sup>* ln <sup>Λ</sup><sup>3</sup>

 ,

. (13)

*u*(*rij*). (14)

 *d***r**1*d***r**2.

*<sup>∂</sup><sup>r</sup> <sup>g</sup>*(*r*)*dr*. (12)

where the mean value is expressed with the pair correlation function by:

<sup>=</sup> *<sup>ρ</sup>*<sup>2</sup> (*<sup>N</sup>* <sup>−</sup> <sup>2</sup>)! *N*!

 *r*<sup>12</sup>

**4.3 Chemical potential and entropy**

introduce a new particle in the system:

From footnote (1), the free energy *F* is written:

so that the chemical potential can be simplified as:

*F*(*V*, *T*, *N*) = −*kBT* ln *QN*(*V*, *T*) = −*kBT*

*μ* = *kBT*

*U*(**r***N*, *λ*) = *λ*

*g*(3)(**r**1,**r**2,**r**3).

above relation:

*∂u*(*r*12) *∂r*<sup>12</sup>

**r** = **r**<sup>1</sup> − **r**2, and simplify the expression of pressure as:

*<sup>p</sup>* <sup>=</sup> *kBT <sup>N</sup>*

*<sup>V</sup>* <sup>−</sup> <sup>2</sup>*<sup>π</sup>* 3 *ρ*2

integrating the partition function with respect to a parameter *λ* to be defined (8).

*μ* = *F*(*V*, *T*, *N*) − *F*(*V*, *T*, *N* − 1) =

<sup>−</sup> ln *ZN*(*V*, *<sup>T</sup>*)

*N* ∑ *j*=2

Secondly, the procedure requires to write the potential energy as a function of the *coupling parameter λ*, under the following form, in order to assess the argument of the logarithm in the

*u*(*r*1*j*) +

Varying from 0 to 1, the coupling parameter *λ* measures the degree of coupling of the particle to which it is assigned (1 in this case) with the rest of the system. In the previous relation, *λ* = 1

$$\mathcal{U}(\mathbf{r}^N, 0) = \sum\_{i}^{N} \sum\_{j>i \ge 2}^{N} \mu(r\_{ij}) = \mathcal{U}(\mathbf{r}^{N-1}).$$

Under these conditions, the configuration integrals for a total coupling (*λ* = 1) and a zero coupling (*λ* = 0) are respectively:

$$Z\_N(V, T, \lambda = 1) = \int \int\_{3N} \exp\left[-\beta \mathcal{U}(\mathbf{r}^N)\right] d\mathbf{r}\_1 d\mathbf{r}\_2...d\mathbf{r}\_N = Z\_N(V, T), \tag{15}$$

$$Z\_N(V, T, \lambda = 0) = \int\_V d\mathbf{r}\_1 \int \int\_{\Im(N-1)} \exp\left[-\beta U(\mathbf{r}^{N-1})\right] d\mathbf{r}\_2...d\mathbf{r}\_N = V Z\_{N-1}(V, T). \tag{16}$$

These expressions are then used to calculate the logarithm of the ratio of configuration integrals in equation (13):

$$\ln \frac{Z\_N(V, T)}{Z\_{N-1}(V, T)} = \ln \frac{Z\_N(V, T, \lambda = 1)}{Z\_N(V, T, \lambda = 0)} + \ln V \tag{17}$$

$$=\ln V + \int\_0^1 \frac{\partial \ln Z\_N}{\partial \lambda} d\lambda. \tag{18}$$

But with the configuration integral *ZN*(*V*, *T*, *λ*), in which potential energy is given by equation (14), we can easily evaluate the partial derivatives *<sup>∂</sup>ZN ∂λ* and *<sup>∂</sup>* ln *ZN ∂λ* . In particular, with the result of the footnote (1), we can write *<sup>∂</sup>* ln *ZN ∂λ* as a function of the pair correlation function as:

$$\frac{\partial \ln Z\_N(V, T, \lambda)}{\partial \lambda} = -\beta \rho^2 \frac{(N - 1)(N - 2)!}{N!} \int \int\_{\mathcal{S}} u(r\_{12}) \left\{ g\_N^{(2)}(\mathbf{r}\_1, \mathbf{r}\_{2'}, \lambda) \right\} d\mathbf{r}\_1 d\mathbf{r}\_2.$$

In addition, if the fluid is homogeneous and isotropic, the above relation simplifies under the following form:

$$\frac{\partial \ln Z\_N(V, T, \lambda)}{\partial \lambda} = -\frac{\beta \rho^2}{N} V \int\_0^\infty u(r) g(r, \lambda) 4\pi r^2 dr$$

that remains only to be substituted in equation (18) for obtaining the logarithm of the ratio of configuration integrals. And by putting the last expression in equation (13), one ultimately arrives to the expression of the chemical potential:

$$
\mu = k\_B T \ln \rho \Lambda^3 + 4 \pi \rho \int\_0^1 \int\_0^\infty u(r) g(r, \lambda) r^2 dr d\lambda. \tag{19}
$$

of Simple Liquids 15

Thermodynamic Perturbation Theory of Simple Liquids 853

Fig. 4. Representation of the hard-sphere potential and its Boltzmann factor.

pressure *pex* comes from the interactions between particles. They are written:

*<sup>ρ</sup>kBT* <sup>=</sup> <sup>4</sup>*<sup>η</sup>* <sup>+</sup> *<sup>η</sup>*2*<sup>B</sup>*

spherical particles on the total volume *V* of the system, that is to say:

*<sup>η</sup>* <sup>=</sup> <sup>1</sup> *V* 4*π* 3 *σ* 2 3

*<sup>ρ</sup>kBT* <sup>4</sup>*<sup>η</sup>* <sup>+</sup> <sup>10</sup>*η*<sup>2</sup> <sup>+</sup> <sup>18</sup>*η*<sup>3</sup> <sup>+</sup> <sup>28</sup>*η*<sup>4</sup> <sup>+</sup> <sup>40</sup>*η*<sup>5</sup> <sup>+</sup> <sup>54</sup>*η*6...

*pex*

and

according to the expansion:

*pex*

The first term of the last equality represents the contribution of the ideal gas, and the excess

where *η* is the *packing fraction* defined by the ratio of the volume actually occupied by the *N*

Note that the first 6 coefficients of the excess pressure *pex* have been calculated analytically and by molecular dynamics (10), with great accuracy. In addition, Carnahan and Starling (11) have shown that the excess pressure of the hard-sphere fluid can be very well predicted by rounding the numerical values of the 6 coefficients towards the nearest integer values,

<sup>3</sup>(*T*) + *<sup>η</sup>*3*<sup>B</sup>*

*<sup>N</sup>* <sup>=</sup> *<sup>π</sup>*

<sup>4</sup>(*T*) + ...

<sup>6</sup> *ρσ*3. (22)

(*k*<sup>2</sup> + 3*k*)*ηk*. (23)

∞ ∑ *k*=1

*pGP <sup>ρ</sup>kBT* <sup>=</sup> 1,

Thus, like the internal energy (Eq. 9) and pressure (Eq. 12), the chemical potential (Eq. 19) is calculated using the pair potential and pair correlation function.

Finally, one writes the entropy *S* in terms of the pair potential and pair correlation function, owing to the expressions of the internal energy (Eq. 9), pressure (Eq. 12) and chemical potential (Eq. 19) (cf. footnote 1):

$$S = \frac{E - F}{T} = \frac{E}{T} - \frac{\mu N}{T} + \frac{pV}{T} \,. \tag{20}$$

It should be noted that the entropy can also be estimated only with the pair correlation function *g*(*r*), without recourse to the pair potential *u*(*r*). The reader interested by this issue should refer to the original articles (9).

#### **4.4 Application to the hard-sphere potential**

In this subsection we determine the equation of state of the hard-sphere system, of which the pair potential being:

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

where *σ* is the hard-sphere diameter. The Boltzmann factor associated with this potential has a significant feature that enable us to express the thermodynamic properties under particularly simple forms. Indeed, the representation of the Boltzmann factor

$$\exp\left[-\beta\mu(r)\right] = \begin{cases} 0 & \text{if} \quad r < \sigma \\ 1 & \text{if} \quad r > \sigma \end{cases}$$

is a step function (Fig. 4) whose derivative with respect to *r* is the Dirac delta function, i. e.:

$$\frac{\partial}{\partial r}\exp\left[-\beta\mu(r)\right] = -\beta\frac{\partial\mu}{\partial r}\exp\left[-\beta\mu(r)\right] = \delta(r-\sigma).$$

In substituting *<sup>∂</sup><sup>u</sup> <sup>∂</sup><sup>r</sup>* , taken from the previous relation, in equation (12) we find the expression of the pressure:

$$p = k\_B T \frac{N}{V} - \frac{2\pi}{3} \rho^2 \int\_0^\infty r^3 \left\{ -\frac{1}{\beta} \frac{\delta(r - \sigma)}{\exp\left[ -\beta \mu(r) \right]} \right\} g(r) dr, \quad \sigma$$

or:

$$p = k\_B T \frac{N}{V} + \frac{2\pi}{3} k\_B T \rho^2 \sigma^3 g(\sigma) \exp\left[\beta u(\sigma)\right]. \tag{21}$$

It is important to recall that, for moderately dense gases, the pressure is usually expressed under the form of the virial expansion

$$\frac{p}{\rho k\_B T} = 1 + \rho B\_2(T) + \rho^2 B\_3(T) + \rho^3 B\_4(T) + \dots = \frac{p^{GP}}{\rho k\_B T} + \frac{p^{ex}}{\rho k\_B T}.$$

Fig. 4. Representation of the hard-sphere potential and its Boltzmann factor.

The first term of the last equality represents the contribution of the ideal gas, and the excess pressure *pex* comes from the interactions between particles. They are written:

$$\frac{p^{GP}}{\rho k\_B T} = 1\_\prime$$

and

14 Thermodynamics book 1

Thus, like the internal energy (Eq. 9) and pressure (Eq. 12), the chemical potential (Eq. 19) is

Finally, one writes the entropy *S* in terms of the pair potential and pair correlation function, owing to the expressions of the internal energy (Eq. 9), pressure (Eq. 12) and chemical

It should be noted that the entropy can also be estimated only with the pair correlation function *g*(*r*), without recourse to the pair potential *u*(*r*). The reader interested by this issue

In this subsection we determine the equation of state of the hard-sphere system, of which the

where *σ* is the hard-sphere diameter. The Boltzmann factor associated with this potential has a significant feature that enable us to express the thermodynamic properties under particularly

> ⎧ ⎨ ⎩

is a step function (Fig. 4) whose derivative with respect to *r* is the Dirac delta function, i. e.:

� ∞ 0 *r* 3 � − 1 *β*

2*π*

*<sup>ρ</sup>kBT* <sup>=</sup> <sup>1</sup> <sup>+</sup> *<sup>ρ</sup>B*2(*T*) + *<sup>ρ</sup>*2*B*3(*T*) + *<sup>ρ</sup>*3*B*4(*T*) + ... <sup>=</sup> *<sup>p</sup>GP*

It is important to recall that, for moderately dense gases, the pressure is usually expressed

∞ if *r* < *σ*

0 if *r* > *σ*,

0 if *r* < *σ*

1 if *r* > *σ*,

*<sup>∂</sup><sup>r</sup>* , taken from the previous relation, in equation (12) we find the expression of

*<sup>∂</sup><sup>r</sup>* exp [−*βu*(*r*)] <sup>=</sup> *<sup>δ</sup>*(*<sup>r</sup>* <sup>−</sup> *<sup>σ</sup>*).

*δ*(*r* − *σ*) exp [−*βu*(*r*)] �

<sup>3</sup> *kBTρ*2*σ*3*g*(*σ*) exp [*βu*(*σ*)] . (21)

*g*(*r*)*dr*,

*<sup>ρ</sup>kBT* <sup>+</sup> *<sup>p</sup>ex*

*<sup>ρ</sup>kBT* .

⎧ ⎨ ⎩ *<sup>T</sup>* <sup>−</sup> *<sup>μ</sup><sup>N</sup>*

*<sup>T</sup>* <sup>+</sup> *pV*

*<sup>T</sup>* . (20)

*<sup>T</sup>* <sup>=</sup> *<sup>E</sup>*

calculated using the pair potential and pair correlation function.

*<sup>S</sup>* <sup>=</sup> *<sup>E</sup>* <sup>−</sup> *<sup>F</sup>*

*u*(*r*) =

exp [−*βu*(*r*)] =

*<sup>∂</sup><sup>r</sup>* exp [−*βu*(*r*)] <sup>=</sup> <sup>−</sup>*<sup>β</sup> <sup>∂</sup><sup>u</sup>*

*<sup>V</sup>* <sup>−</sup> <sup>2</sup>*<sup>π</sup>* 3 *ρ*2

*<sup>p</sup>* <sup>=</sup> *kBT <sup>N</sup>*

*<sup>V</sup>* <sup>+</sup>

simple forms. Indeed, the representation of the Boltzmann factor

potential (Eq. 19) (cf. footnote 1):

should refer to the original articles (9).

pair potential being:

In substituting *<sup>∂</sup><sup>u</sup>*

the pressure:

or:

**4.4 Application to the hard-sphere potential**

*∂*

*<sup>p</sup>* <sup>=</sup> *kBT <sup>N</sup>*

under the form of the virial expansion

*p*

$$\frac{p^{\varepsilon\chi}}{\rho k\_B T} = 4\eta + \eta^2 B\_3'(T) + \eta^3 B\_4'(T) + \dots$$

where *η* is the *packing fraction* defined by the ratio of the volume actually occupied by the *N* spherical particles on the total volume *V* of the system, that is to say:

$$
\eta = \frac{1}{V} \frac{4\pi}{3} \left(\frac{\sigma}{2}\right)^3\\N = \frac{\pi}{6} \rho \sigma^3. \tag{22}
$$

Note that the first 6 coefficients of the excess pressure *pex* have been calculated analytically and by molecular dynamics (10), with great accuracy. In addition, Carnahan and Starling (11) have shown that the excess pressure of the hard-sphere fluid can be very well predicted by rounding the numerical values of the 6 coefficients towards the nearest integer values, according to the expansion:

$$\frac{p^{ex}}{\rho k\_B T} \simeq 4\eta + 10\eta^2 + 18\eta^3 + 28\eta^4 + 40\eta^5 + 54\eta^6 \dots \simeq \sum\_{k=1}^{\infty} (k^2 + 3k)\eta^k. \tag{23}$$

of Simple Liquids 17

Thermodynamic Perturbation Theory of Simple Liquids 855

Like the pressure, this expansion is written as a rational function by combining the geometric

(*<sup>k</sup>* <sup>+</sup> <sup>3</sup>)*η<sup>k</sup>* <sup>=</sup> <sup>4</sup>*<sup>η</sup>* <sup>−</sup> <sup>3</sup>*η*<sup>2</sup>

*NkBT* <sup>=</sup> ln *<sup>ρ</sup>*Λ<sup>3</sup> <sup>−</sup> <sup>1</sup> <sup>+</sup>

ln *<sup>ρ</sup>*Λ<sup>3</sup> <sup>−</sup> <sup>5</sup>

= −*NkB*

5

Finally, combining equations (25) and (24), with the help of equation (20), one reaches the

*<sup>ρ</sup>kBT* <sup>=</sup> ln *<sup>ρ</sup>*Λ<sup>3</sup> <sup>−</sup> <sup>1</sup> <sup>+</sup>

∞ ∑ *k*=1

(1 − *η*)

*kη<sup>k</sup>* + ∞ ∑ *k*=1 3*η<sup>k</sup>* ,

<sup>2</sup> +

3*η* (1 − *η*) ,

(*k* + 3)*η<sup>k</sup>* =

(*<sup>k</sup>* <sup>+</sup> <sup>3</sup>)*η<sup>k</sup>* <sup>=</sup> *<sup>η</sup>*

<sup>2</sup> <sup>−</sup> <sup>4</sup>*<sup>η</sup>* <sup>−</sup> <sup>3</sup>*η*<sup>2</sup> (1 − *η*)

2 ,

<sup>4</sup>*<sup>η</sup>* <sup>−</sup> <sup>3</sup>*η*<sup>2</sup> (1 − *η*)

2 ,

<sup>1</sup> <sup>+</sup> <sup>5</sup>*<sup>η</sup>* <sup>−</sup> <sup>6</sup>*η*<sup>2</sup> <sup>+</sup> <sup>2</sup>*η*<sup>3</sup> (1 − *η*)

<sup>2</sup> . (26)

<sup>3</sup> . (27)

*V* −

(1 − *η*)

2 ,

 *∂Fex ∂T*

<sup>4</sup>*<sup>η</sup>* <sup>−</sup> <sup>3</sup>*η*<sup>2</sup> (1 − *η*)

> *V* ,

<sup>2</sup> . (25)

*<sup>k</sup>*=<sup>1</sup> *<sup>η</sup><sup>k</sup>* with its first derivative6. The expression of the excess free energy is:

*<sup>η</sup>* , *<sup>F</sup>ex* is then reduced to the series expansion:

∞ ∑ *k*=1

(*k* + 3)*ηk*.

*<sup>d</sup><sup>η</sup>* <sup>=</sup> <sup>−</sup>*<sup>V</sup>*

*NkBT* <sup>=</sup> <sup>4</sup>*<sup>η</sup>* <sup>+</sup> <sup>5</sup>*η*<sup>2</sup> <sup>+</sup> <sup>6</sup>*η*<sup>3</sup> <sup>+</sup> <sup>7</sup>*η*<sup>4</sup> <sup>+</sup> <sup>8</sup>*η*<sup>5</sup> <sup>+</sup> <sup>9</sup>*η*6... <sup>=</sup>

∞ ∑ *k*=1

where *SGP* is the entropy of the ideal gas given by the Sackur-Tetrode equation:

*<sup>S</sup>GP* <sup>=</sup> <sup>−</sup>*NkB*

 *∂Fex ∂T*

 *V*

<sup>=</sup> <sup>−</sup> ln *<sup>ρ</sup>*Λ<sup>3</sup> <sup>+</sup>

*Fex*

Now, the entropy is obtained using the same method of calculation, by deriving the free

and the free energy of the hard-sphere fluid reduces to the following form:

*NkBT* <sup>+</sup>

 *∂F ∂T V* = − *∂FGP ∂T*

*Fex NkBT* <sup>=</sup>

But, with equation (22) that gives *dV*

*Fex*

*F NkBT* <sup>=</sup> *<sup>F</sup>GP*

*S* = −

and where the excess entropy *Sex* arises from the relation:

*<sup>S</sup>ex* <sup>=</sup> <sup>−</sup>

*S NkB*

*NkBT* <sup>+</sup> *<sup>p</sup>*

∞ ∑ *k*=1

∞ ∑ *k*=1

hence the expression of the entropy of the hard-sphere fluid:

chemical potential of the hard-sphere fluid that reads:

*μ kBT* <sup>=</sup> *<sup>F</sup>*

<sup>6</sup> Indeed, the identity:

is yet written:

energy with respect to temperature:

series ∑<sup>∞</sup>

In combining the first and second derivatives of the geometric series ∑<sup>∞</sup> *<sup>k</sup>*=<sup>1</sup> *<sup>η</sup>k*, it is found that equation (23) can be transformed into a rational fraction5 enabling the deduction of the excess pressure in the form:

$$\frac{p^{\varepsilon\chi}}{\rho k\_B T} \simeq \sum\_{k=1}^{\infty} (k^2 + 3k)\eta^k = \frac{4\eta - 2\eta^2}{(1-\eta)^3}.$$

Consequently, the equation of state of the hard-sphere fluid is written with excellent precision as:

$$\frac{p}{\rho k\_B T} = \frac{1 + \eta + \eta^2 - \eta^3}{\left(1 - \eta\right)^3}.\tag{24}$$

It is also possible to calculate the internal energy of the hard-sphere fluid by substituting *u*(*r*) in equation (9). Given that *u*(*r*) is zero when *r* > *σ* and *g*(*r*) is zero when *r* < *σ*, it follows that the integral is always zero, and that the internal energy of the hard-sphere fluid is equal to that of the ideal gas *E* = <sup>3</sup> <sup>2</sup> *NkBT*.

As for the free energy *F*, it is determined by integrating the pressure over volume with the equation:

$$p = -\left(\frac{\partial F}{\partial V}\right)\_T = -\left(\frac{\partial F^{GP}}{\partial V}\right)\_T - \left(\frac{\partial F^{\epsilon \chi}}{\partial V}\right)\_T$$

where *FGP* is the free energy of ideal gas (cf. footnote 1, with *ZN*(*V*, *T*) = *VN*):

$$F^{GP} = Nk\_B T \left(\ln \rho \Lambda^3 - 1\right),$$

and *Fex* the excess free energy, calculated by integrating equation (23) as follows:

$$F^{\rm ex} = -\int p^{\rm ex}dV = -\int \frac{Nk\_BT}{V} \left( 4\eta + 10\eta^2 + 18\eta^3 + 28\eta^4 + 40\eta^5 + 54\eta^6... \right) \frac{dV}{d\eta} d\eta.$$

<sup>5</sup> To obtain the rational fraction, one must decompose the sum as:

$$\sum\_{k=1}^{\infty} (k^2 + 3k)\eta^k = \sum\_{k=1}^{\infty} (k^2 - k)\eta^k + \sum\_{k=1}^{\infty} 4k\eta^k.$$

and combine the geometric series ∑<sup>∞</sup> *<sup>k</sup>*=<sup>1</sup> *<sup>η</sup><sup>k</sup>* with its first and second derivatives:

$$\begin{aligned} \sum\_{k=1}^{\infty} \eta^k &= \eta + \eta^2 + \eta^3 + \dots = \frac{\eta}{1-\eta}, \\ \sum\_{k=1}^{\infty} k \eta^{k-1} &= \frac{1}{\left(1-\eta\right)^2}, \\ \sum\_{k=1}^{\infty} k(k-1)\eta^{k-2} &= \frac{2}{\left(1-\eta\right)^3}. \end{aligned}$$

and

$$\text{to see appear the relation:}$$

$$\sum\_{k=1}^{\infty} (k^2 + 3k)\eta^k = \frac{2\eta^2}{\left(1 - \eta\right)^3} + \frac{4\eta}{\left(1 - \eta\right)^2}.$$

16 Thermodynamics book 1

equation (23) can be transformed into a rational fraction5 enabling the deduction of the excess

Consequently, the equation of state of the hard-sphere fluid is written with excellent precision

*<sup>ρ</sup>kBT* <sup>=</sup> <sup>1</sup> <sup>+</sup> *<sup>η</sup>* <sup>+</sup> *<sup>η</sup>*<sup>2</sup> <sup>−</sup> *<sup>η</sup>*<sup>3</sup> (1 − *η*)

It is also possible to calculate the internal energy of the hard-sphere fluid by substituting *u*(*r*) in equation (9). Given that *u*(*r*) is zero when *r* > *σ* and *g*(*r*) is zero when *r* < *σ*, it follows that the integral is always zero, and that the internal energy of the hard-sphere fluid is equal

As for the free energy *F*, it is determined by integrating the pressure over volume with the

∞ ∑ *k*=1

*<sup>k</sup>ηk*−<sup>1</sup> <sup>=</sup> <sup>1</sup>

*<sup>k</sup>*(*<sup>k</sup>* <sup>−</sup> <sup>1</sup>)*ηk*−<sup>2</sup> <sup>=</sup> <sup>2</sup>

(*k*<sup>2</sup> <sup>+</sup> <sup>3</sup>*k*)*η<sup>k</sup>* <sup>=</sup> <sup>2</sup>*η*<sup>2</sup>

∞ ∑ *k*=1

∞ ∑ *k*=1

∞ ∑ *k*=1

∞ ∑ *k*=1 ln *<sup>ρ</sup>*Λ<sup>3</sup> <sup>−</sup> <sup>1</sup>

(*k*<sup>2</sup> <sup>−</sup> *<sup>k</sup>*)*η<sup>k</sup>* <sup>+</sup>

*<sup>k</sup>*=<sup>1</sup> *<sup>η</sup><sup>k</sup>* with its first and second derivatives:

(1 − *η*) 2 ,

(1 − *η*) 3 ,

<sup>3</sup> +

4*η* (1 − *η*) 2 .

(1 − *η*)

*<sup>η</sup><sup>k</sup>* <sup>=</sup> *<sup>η</sup>* <sup>+</sup> *<sup>η</sup>*<sup>2</sup> <sup>+</sup> *<sup>η</sup>*<sup>3</sup> <sup>+</sup> ... <sup>=</sup> *<sup>η</sup>*

*T* −  *∂Fex ∂V*

 ,

4*η* + 10*η*<sup>2</sup> + 18*η*<sup>3</sup> + 28*η*<sup>4</sup> + 40*η*<sup>5</sup> + 54*η*6...

∞ ∑ *k*=1 4*kη<sup>k</sup>* ,

<sup>1</sup> <sup>−</sup> *<sup>η</sup>* ,

 *T* ,

(*k*<sup>2</sup> <sup>+</sup> <sup>3</sup>*k*)*η<sup>k</sup>* <sup>=</sup> <sup>4</sup>*<sup>η</sup>* <sup>−</sup> <sup>2</sup>*η*<sup>2</sup>

(1 − *η*)

3 .

<sup>3</sup> . (24)

*<sup>k</sup>*=<sup>1</sup> *<sup>η</sup>k*, it is found that

 *dV dη dη*.

In combining the first and second derivatives of the geometric series ∑<sup>∞</sup>

∞ ∑ *k*=1

*p*

*pex <sup>ρ</sup>kBT* �

<sup>2</sup> *NkBT*.

 *∂F ∂V T* = − *∂FGP ∂V*

 *NkBT V*

<sup>5</sup> To obtain the rational fraction, one must decompose the sum as: ∞ ∑ *k*=1

where *FGP* is the free energy of ideal gas (cf. footnote 1, with *ZN*(*V*, *T*) = *VN*):

*FGP* = *NkBT*

and *Fex* the excess free energy, calculated by integrating equation (23) as follows:

(*k*<sup>2</sup> + 3*k*)*η<sup>k</sup>* =

*p* = −

pressure in the form:

to that of the ideal gas *E* = <sup>3</sup>

as:

equation:

*<sup>F</sup>ex* <sup>=</sup> <sup>−</sup>

*<sup>p</sup>exdV* <sup>=</sup> <sup>−</sup>

and combine the geometric series ∑<sup>∞</sup>

to see appear the relation:

and

But, with equation (22) that gives *dV <sup>d</sup><sup>η</sup>* <sup>=</sup> <sup>−</sup>*<sup>V</sup> <sup>η</sup>* , *<sup>F</sup>ex* is then reduced to the series expansion:

$$\frac{F^{ex}}{Nk\_BT} = 4\eta + 5\eta^2 + 6\eta^3 + 7\eta^4 + 8\eta^5 + 9\eta^6 \dots = \sum\_{k=1}^{\infty} (k+3)\eta^k X$$

Like the pressure, this expansion is written as a rational function by combining the geometric series ∑<sup>∞</sup> *<sup>k</sup>*=<sup>1</sup> *<sup>η</sup><sup>k</sup>* with its first derivative6. The expression of the excess free energy is:

$$\frac{F^{\mathrm{ex}}}{Nk\_B T} = \sum\_{k=1}^{\infty} (k+3)\eta^k = \frac{4\eta - 3\eta^2}{\left(1 - \eta\right)^2}.$$

and the free energy of the hard-sphere fluid reduces to the following form:

$$\frac{F}{Nk\_BT} = \frac{F^{GP}}{Nk\_BT} + \frac{F^{ex}}{Nk\_BT} = \ln\rho\Lambda^3 - 1 + \frac{4\eta - 3\eta^2}{\left(1 - \eta\right)^2}.\tag{25}$$

Now, the entropy is obtained using the same method of calculation, by deriving the free energy with respect to temperature:

$$S = -\left(\frac{\partial F}{\partial T}\right)\_V = -\left(\frac{\partial F^{GP}}{\partial T}\right)\_V - \left(\frac{\partial F^{\text{ex}}}{\partial T}\right)\_V$$

where *SGP* is the entropy of the ideal gas given by the Sackur-Tetrode equation:

$$S^{GP} = -Nk\_B \left(\ln \rho \Lambda^3 - \frac{5}{2}\right) \Lambda$$

and where the excess entropy *Sex* arises from the relation:

$$S^{\varepsilon \chi} = -\left(\frac{\partial F^{\varepsilon \chi}}{\partial T}\right)\_V = -Nk\_B \frac{4\eta - 3\eta^2}{\left(1 - \eta\right)^2}.$$

hence the expression of the entropy of the hard-sphere fluid:

$$\frac{S}{Nk\_B} = -\ln \rho \Lambda^3 + \frac{5}{2} - \frac{4\eta - 3\eta^2}{\left(1 - \eta\right)^2}. \tag{26}$$

,

Finally, combining equations (25) and (24), with the help of equation (20), one reaches the chemical potential of the hard-sphere fluid that reads:

$$\frac{\mu}{k\_B T} = \frac{F}{Nk\_B T} + \frac{p}{\rho k\_B T} = \ln \rho \Lambda^3 - 1 + \frac{1 + 5\eta - 6\eta^2 + 2\eta^3}{\left(1 - \eta\right)^3}. \tag{27}$$

<sup>6</sup> Indeed, the identity:

$$
\sum\_{k=1}^{\infty} (k+3)\eta^k = \sum\_{k=1}^{\infty} k\eta^k + \sum\_{k=1}^{\infty} 3\eta^k,
$$

$$
\sum\_{k=1}^{\infty} (k+3)\eta^k = \frac{\eta}{(1-\eta)^2} + \frac{3\eta}{(1-\eta)^3}
$$

of Simple Liquids 19

Thermodynamic Perturbation Theory of Simple Liquids 857

Such a qualitative success emphasizes the role played by the repulsive part of the pair potential to describe the structure factor of liquids, while the long-ranged attractive contribution has a minor role. It can be said for simplicity that the repulsive contribution of the potential determines the structure of liquids (stacking of atoms and steric effects) and

It is important to remember that the thermodynamic properties of the hard-sphere fluid (Eqs. 24, 25, 26, 27) and the structure factor *SHS*(*q*) can be calculated with great accuracy. That suggests replacing the repulsive part of potential in real systems by the hard-sphere potential that becomes the *reference system*, and precict the structural and thermodynamic properties of real systems with those of the hard-sphere fluid, after making the necessary adaptations. To perform these adaptations, the attractive contribution of potential should be treated as a

The rest of this subsection is devoted to a summary of thermodynamic perturbation methods7. It should be noted, from the outset, that the calculation of thermodynamic properties with the thermodynamic perturbation methods requires knowledge of the pair correlation function

In perturbation theory proposed by Zwanzig (16), it is assumed that the total potential energy *U*(**r***N*) of the system can be divided into two parts. The first part, *U*0(**r***N*), is the energy of the unperturbed system considered as reference system and the second part, *U*1(**r***N*), is the energy of the perturbation which is much smaller that *U*0(**r***N*). More precisely, it is posed that

*U*(**r***N*) = *U*0(**r***N*) + *λU*1(**r***N*)

in order to vary continuously the potential energy from *U*0(**r***N*) to *U*(**r***N*), by changing *λ* from

By replacing the potential energy *U*(**r***N*) in the expression of the configuration integral (cf.

 3*N* exp

 *∂*2*F ∂λ*<sup>2</sup> 

<sup>−</sup>*βλU*1(**r***N*)

 0

<sup>−</sup>*βλU*1(**r***N*)

+ ... (28)

<sup>−</sup>*βU*0(**r***N*)

<sup>−</sup>*βλU*1(**r***N*)

, (29)

 *d***r***<sup>N</sup>*

> , so

exp

<sup>3</sup>*<sup>N</sup>* exp [−*βU*0(**r***N*)] *<sup>d</sup>***r***<sup>N</sup>* .

*<sup>N</sup>* (*V*, *T*) of the reference system, and

the attractive contribution is responsible for their cohesion.

*gHS*(*r*) of the hard-sphere system and not that of the real system.

the potential energy depend on the coupling parameter *λ* by the relation:

0 to 1, and that the free energy *F* of the system is expanded in Taylor series as:

 *∂F ∂λ* + *λ*2 2

*F* = *F*<sup>0</sup> + *λ*

 *<sup>d</sup>***r***<sup>N</sup>* <sup>×</sup>

the remaining term can be regarded as the average value of the quantity exp

*<sup>N</sup>* (*V*, *T*)

 exp

where �...�<sup>0</sup> refers to the statistical average in the canonical ensemble of the reference system. After the substitution of the configuration integral (Eq. 29) in the expression of the free energy

<sup>7</sup> The interested reader will find all useful adjuncts in the books either by J. P. Hansen and I. R. McDonald

<sup>−</sup>*βU*0(**r***N*)

The first integral represents the configuration integral *Z*(0)

that the previous relation can be put under the general form:

*ZN*(*V*, *<sup>T</sup>*) = *<sup>Z</sup>*(0)

perturbation to the reference system.

**5.1 Zwanzig method**

footnote 1), one gets:

or by D. A. McQuarrie.

3*N*

exp

*ZN*(*V*, *<sup>T</sup>*) =

Since they result from equation (23), the expressions of thermodynamic properties (*p*, *F*, *S* and *μ*) of the hard-sphere fluid make up a homogeneous group of relations related to the Carnahan and Starling equation of state. But other expressions of thermodynamic properties can also be determined using the pressure equation of state (Eq. 12) and the compressibility equation of state, which will not be discussed here. Unlike the Carnahan and Starling equation of state, these two equations of state require knowledge of the pair correlation function of hard spheres, *gHS*(*r*). The latter is not available in analytical form. The interested reader will find the Fortran program aimed at doing its calculation, in the book by McQuarrie (12), page 600. It should be mentioned that the thermodynamic properties (*p*, *F*, *S* and *μ*), obtained with the equations of state of pressure and compressibility, have analytical forms similar to those from the Carnahan and Starling equation of state, and they provide results whose differences are indistinguishable to low densities.

### **5. Thermodynamic perturbation theory**

All theoretical and experimental studies have shown that the structure factor *S*(*q*) of simple liquids resembles that of the hard-sphere fluid. For proof, just look at the experimental structure factor of liquid sodium (13) at 373 K, in comparison with the structure factor of hard-sphere fluid (14) for a value of the packing fraction *η* of 0.45. We can see that the agreement is not bad, although there is a slight shift of the oscillations and ratios of peak heights significantly different. Besides, numerical calculations showed that the structure factor obtained with the Lennard-Jones potential describes the structure of simple fluids (15) and looks like the structure factor of hard-sphere fluid whose diameter is chosen correctly.

Fig. 5. Experimental structure factor of liquid sodium at 373 K (points), and hard-sphere structure factor (solid curve), with *η* = 0, 45.

Such a qualitative success emphasizes the role played by the repulsive part of the pair potential to describe the structure factor of liquids, while the long-ranged attractive contribution has a minor role. It can be said for simplicity that the repulsive contribution of the potential determines the structure of liquids (stacking of atoms and steric effects) and the attractive contribution is responsible for their cohesion.

It is important to remember that the thermodynamic properties of the hard-sphere fluid (Eqs. 24, 25, 26, 27) and the structure factor *SHS*(*q*) can be calculated with great accuracy. That suggests replacing the repulsive part of potential in real systems by the hard-sphere potential that becomes the *reference system*, and precict the structural and thermodynamic properties of real systems with those of the hard-sphere fluid, after making the necessary adaptations. To perform these adaptations, the attractive contribution of potential should be treated as a perturbation to the reference system.

The rest of this subsection is devoted to a summary of thermodynamic perturbation methods7. It should be noted, from the outset, that the calculation of thermodynamic properties with the thermodynamic perturbation methods requires knowledge of the pair correlation function *gHS*(*r*) of the hard-sphere system and not that of the real system.

#### **5.1 Zwanzig method**

18 Thermodynamics book 1

Since they result from equation (23), the expressions of thermodynamic properties (*p*, *F*, *S* and *μ*) of the hard-sphere fluid make up a homogeneous group of relations related to the Carnahan and Starling equation of state. But other expressions of thermodynamic properties can also be determined using the pressure equation of state (Eq. 12) and the compressibility equation of state, which will not be discussed here. Unlike the Carnahan and Starling equation of state, these two equations of state require knowledge of the pair correlation function of hard spheres, *gHS*(*r*). The latter is not available in analytical form. The interested reader will find the Fortran program aimed at doing its calculation, in the book by McQuarrie (12), page 600. It should be mentioned that the thermodynamic properties (*p*, *F*, *S* and *μ*), obtained with the equations of state of pressure and compressibility, have analytical forms similar to those from the Carnahan and Starling equation of state, and they provide results whose differences are

All theoretical and experimental studies have shown that the structure factor *S*(*q*) of simple liquids resembles that of the hard-sphere fluid. For proof, just look at the experimental structure factor of liquid sodium (13) at 373 K, in comparison with the structure factor of hard-sphere fluid (14) for a value of the packing fraction *η* of 0.45. We can see that the agreement is not bad, although there is a slight shift of the oscillations and ratios of peak heights significantly different. Besides, numerical calculations showed that the structure factor obtained with the Lennard-Jones potential describes the structure of simple fluids (15) and looks like the structure factor of hard-sphere fluid whose diameter is chosen correctly.

Fig. 5. Experimental structure factor of liquid sodium at 373 K (points), and hard-sphere

indistinguishable to low densities.

**5. Thermodynamic perturbation theory**

structure factor (solid curve), with *η* = 0, 45.

In perturbation theory proposed by Zwanzig (16), it is assumed that the total potential energy *U*(**r***N*) of the system can be divided into two parts. The first part, *U*0(**r***N*), is the energy of the unperturbed system considered as reference system and the second part, *U*1(**r***N*), is the energy of the perturbation which is much smaller that *U*0(**r***N*). More precisely, it is posed that the potential energy depend on the coupling parameter *λ* by the relation:

$$\mathcal{U}(\mathbf{r}^N) = \mathcal{U}\_0(\mathbf{r}^N) + \lambda \mathcal{U}\_1(\mathbf{r}^N)$$

in order to vary continuously the potential energy from *U*0(**r***N*) to *U*(**r***N*), by changing *λ* from 0 to 1, and that the free energy *F* of the system is expanded in Taylor series as:

$$F = F\_0 + \lambda \left(\frac{\partial F}{\partial \lambda}\right) + \frac{\lambda^2}{2} \left(\frac{\partial^2 F}{\partial \lambda^2}\right) + \dots \tag{28}$$

By replacing the potential energy *U*(**r***N*) in the expression of the configuration integral (cf. footnote 1), one gets:

$$Z\_N(V,T) = \int \int\_{3N} \exp\left[-\beta l l\_0(\mathbf{r}^N)\right] d\mathbf{r}^N \times \frac{\int \int\_{3N} \left\{ \exp\left[-\beta \lambda l l\_1(\mathbf{r}^N)\right] \right\} \exp\left[-\beta l l\_0(\mathbf{r}^N)\right] d\mathbf{r}^N}{\int \int\_{3N} \exp\left[-\beta l l\_0(\mathbf{r}^N)\right] d\mathbf{r}^N}.$$

The first integral represents the configuration integral *Z*(0) *<sup>N</sup>* (*V*, *T*) of the reference system, and the remaining term can be regarded as the average value of the quantity exp <sup>−</sup>*βλU*1(**r***N*) , so that the previous relation can be put under the general form:

$$Z\_N(V, T) = Z\_N^{(0)}(V, T) \left\langle \exp\left[-\beta\lambda U\_1(\mathbf{r}^N)\right] \right\rangle\_{0'} \tag{29}$$

where �...�<sup>0</sup> refers to the statistical average in the canonical ensemble of the reference system. After the substitution of the configuration integral (Eq. 29) in the expression of the free energy

<sup>7</sup> The interested reader will find all useful adjuncts in the books either by J. P. Hansen and I. R. McDonald or by D. A. McQuarrie.

of Simple Liquids 21

Thermodynamic Perturbation Theory of Simple Liquids 859

As a first application of thermodynamic perturbation method, we search the phenomenological van der Waals equation of state. In view of this, consider equation (37) at zero order in *β*. The simplest assumption to determine *c*<sup>1</sup> is to admit that the total

> *i* ∑ *j*>*i*

> > exp

Therefore, the free energy of the perturbation to zero order in *β* is given by equation (33), that

∑*<sup>i</sup>* ∑*j*>*<sup>i</sup> u*1(*rij*)

To simplify the above relation, we proceed as for calculating the internal energy of liquids (Eq. 8) by revealing the pair correlation function of the reference system in the numerator. If we assume that the sum of pair potentials is composed of equivalent terms equal to

<sup>2</sup> *u*1(*r*12), the expression of *c*<sup>1</sup> is simplified as:

The integral in between the braces is then expressed as a function of the pair correlation

Yet, to find the equation of van der Waals we have to choose the hard-sphere system of diameter *σ*, as reference system, and suppose that the perturbation is a long-range potential, weakly attractive, the form of which is not useful to specify (Fig. 6a). Since one was unaware of the existence of the pair correlation function when the model was developed by van der Waals, it is reasonable to estimate *g*0(*r*) by a function equal to zero within the particle, and to one at the outside. According to van der Waals, suppose further that the available volume per

With these simplifications in mind, the configuration integral and free energy of the reference

<sup>8</sup> The parameter *b* introduced by van der Waals is the covolume. Its expression comes from the fact that

*N*

*<sup>d</sup>***r***<sup>N</sup>* <sup>=</sup> (*<sup>V</sup>* <sup>−</sup> *Nb*)

= −*NkBT*

*N* ,

ln (*<sup>V</sup>* <sup>−</sup> *Nb*)

*<sup>N</sup>* <sup>−</sup> 3 ln <sup>Λ</sup> <sup>+</sup> <sup>1</sup>

<sup>3</sup>*πσ*<sup>3</sup> must be assigned to each particle (Fig.

 .

3(*N*−2)

exp

*u*0(*rij*) + ∑

<sup>3</sup>*<sup>N</sup>* exp [−*βU*0(**r***N*)] *<sup>d</sup>***r***<sup>N</sup>* .

*i* ∑ *j*>*i*

<sup>−</sup>*βU*0(**r***N*)

<sup>−</sup>*βU*0(**r***N*)

*u*1(*rij*).

 *d***r***<sup>N</sup>*

*u*1(*r*)*g*0(*r*)*d***r**. (38)

*d***r**3...*d***r***<sup>N</sup>*

*d***r**1*d***r**2.

potential energy may be decomposed into a sum of pair potentials in the form:

*<sup>U</sup>*(**r***N*) = *<sup>U</sup>*0(**r***N*) + *<sup>U</sup>*1(**r***N*) = ∑

 3*N* 

<sup>2</sup> *<sup>u</sup>*1(*r*12)

*<sup>c</sup>*<sup>1</sup> <sup>=</sup> *<sup>ρ</sup>*<sup>2</sup> 2 *d***R** 

<sup>3</sup>*πσ*<sup>3</sup> and the unoccupied volume is (*<sup>V</sup>* <sup>−</sup> *Nb*).

<sup>−</sup>*βU*0(**r***N*)

(*V* − *Nb*)

*N*!Λ3*<sup>N</sup>*

*c*<sup>1</sup> = �*U*1�<sup>0</sup> =

*N*(*N* − 1)

function (cf. footnote 2), and *c*<sup>1</sup> reduces to:

system are respectively (cf. footnote 1):

*<sup>Z</sup>*0(*V*, *<sup>T</sup>*) =

and *<sup>F</sup>*<sup>0</sup> <sup>=</sup> <sup>−</sup>*kBT* ln

3*N*

if two particles are in contact, half of the excluded volume <sup>4</sup>

exp

**5.2 Van der Waals equation**

<sup>∑</sup>*<sup>i</sup>* <sup>∑</sup>*j*>*<sup>i</sup> <sup>u</sup>*1(*rij*) = *<sup>N</sup>*(*N*−1)

*Z*0(*V*, *T*)

*<sup>c</sup>*<sup>1</sup> <sup>=</sup> <sup>1</sup>

particle8 is *b* = <sup>2</sup>

6b).

is to say:

(cf. footnote 1), this one reads:

$$-\beta F = \ln \frac{Z\_N^{(0)}(V, T)}{N! \Lambda^{3N}} + \ln \left\langle \exp \left[ -\beta \lambda U\_1(\mathbf{r}^N) \right] \right\rangle\_0. \tag{30}$$

The first term on the RHS stands for the free energy of the reference system, denoted (−*βF*0), and the second term represents the free energy of the perturbation:

$$-\beta F\_1 = \ln \left\langle \exp \left[ -\beta \lambda U\_1(\mathbf{r}^N) \right] \right\rangle\_0. \tag{31}$$

Since the perturbation *<sup>U</sup>*1(**r***N*) is small, exp (−*βλU*1) can be expanded in series, so that the statistical average exp <sup>−</sup>*βλU*1(**r***N*) <sup>0</sup>, calculated on the reference system, is expressed as:

$$\left\langle \exp \left[ -\beta \lambda \mathcal{U}\_1(\mathbf{r}^N) \right] \right\rangle\_0 = 1 - \beta \lambda \left\langle \mathcal{U}\_1 \right\rangle\_0 + \frac{1}{2!} \beta^2 \lambda^2 \left\langle \mathcal{U}\_1^2 \right\rangle\_0 - \frac{1}{3!} \beta^3 \lambda^3 \left\langle \mathcal{U}\_1^3 \right\rangle\_0 + \dots \tag{32}$$

Incidentally, we may note that the coefficients of *β* in the preceding expansion represent statistical moments in the strict sense. Given the shape of equation (32), it is still possible to write equation (31) by expanding ln exp <sup>−</sup>*βλU*1(**r***N*) <sup>0</sup> in Taylor series. After simplifications, equation (31) reduces to:

$$\begin{split} \left\langle \ln \left\langle \exp \left[ -\beta \lambda \mathcal{U}\_{1} (\mathbf{r}^{N}) \right] \right\rangle\_{0} = -\beta \lambda \left\langle \mathcal{U}\_{1} \right\rangle\_{0} + \frac{1}{2!} \beta^{2} \lambda^{2} \left[ \left\langle \mathcal{U}\_{1}^{2} \right\rangle\_{0} - \left\langle \mathcal{U}\_{1} \right\rangle\_{0}^{2} \right] \\ - \beta^{3} \lambda^{3} \left[ \frac{1}{3!} \left\langle \mathcal{U}\_{1}^{3} \right\rangle\_{0} - \frac{1}{2} \left\langle \mathcal{U}\_{1} \right\rangle\_{0} \left\langle \mathcal{U}\_{1}^{2} \right\rangle\_{0} + \frac{1}{3} \left\langle \mathcal{U}\_{1} \right\rangle\_{0}^{3} \right] + \beta^{4} \lambda^{4} \left[ \dots \right] - \dots \end{split}$$

Now if we set:

*c*<sup>1</sup> = �*U*1�<sup>0</sup> , (33)

$$\mathcal{L}\_2 = \frac{1}{2!} \left[ \left< \mathcal{U}\_1^2 \right>\_0 - \left< \mathcal{U}\_1 \right>\_0^2 \right],\tag{34}$$

$$\mathcal{L}\_3 = \frac{1}{3!} \left[ \left< \mathcal{U}\_1^3 \right>\_0 - 3 \left< \mathcal{U}\_1 \right>\_0 \left< \mathcal{U}\_1^2 \right>\_0 + 2 \left< \mathcal{U}\_1 \right>\_0^3 \right], \quad \text{etc.} \tag{35}$$

we find that:

$$\begin{aligned} \left\langle \ln \left\langle \exp \left[ -\beta \lambda U\_1(\mathbf{r}^N) \right] \right\rangle \right\rangle\_0 &= -\lambda \beta \mathfrak{c}\_1 + \lambda^2 \beta^2 \mathfrak{c}\_2 - \lambda^3 \beta^3 \mathfrak{c}\_3 + \dots \\ \left\langle \right\rangle\_{\mathfrak{c}\_1} &= \lambda \mathfrak{c}\_1 + \lambda^2 \mathfrak{c}\_2 - \lambda^3 \beta^3 \mathfrak{c}\_3 + \dots \end{aligned}$$

The contribution of the perturbation (Eq. 31) to the free energy is then written in the compact form:

$$-\beta F\_1 = \ln \left\langle \exp \left[ -\beta \lambda U\_1(\mathbf{r}^N) \right] \right\rangle\_0 = -\lambda \beta \sum\_{n=1}^{\infty} c\_n (-\lambda \beta)^{n-1},\tag{36}$$

and the expression of the free energy *F* of the real system is found by substituting equation (36) into equation (30), as follows:

$$F = F\_0 + F\_1 = F\_0 + \lambda c\_1 - \lambda^2 \beta c\_2 + \lambda^3 \beta^2 c\_3 + \dots \tag{37}$$

where the free energy of the real system is obtained by putting *λ* = 1. This expression of the free energy of liquids in power series expansion of *β* corresponds to the *high temperature approximation*.

#### **5.2 Van der Waals equation**

20 Thermodynamics book 1

The first term on the RHS stands for the free energy of the reference system, denoted (−*βF*0),

Since the perturbation *<sup>U</sup>*1(**r***N*) is small, exp (−*βλU*1) can be expanded in series, so that the

Incidentally, we may note that the coefficients of *β* in the preceding expansion represent statistical moments in the strict sense. Given the shape of equation (32), it is still

> 1 2! *<sup>β</sup>*2*λ*<sup>2</sup>

− 3 �*U*1�<sup>0</sup>

<sup>−</sup>*βλU*1(**r***N*)

 *U*2 1 0

The contribution of the perturbation (Eq. 31) to the free energy is then written in the compact

and the expression of the free energy *F* of the real system is found by substituting equation

where the free energy of the real system is obtained by putting *λ* = 1. This expression of the free energy of liquids in power series expansion of *β* corresponds to the *high temperature*

1 2! *<sup>β</sup>*2*λ*<sup>2</sup>

exp

 *U*2 1 0 − �*U*1� 2 0 

 exp   exp 

<sup>−</sup>*βλU*1(**r***N*)

<sup>−</sup>*βλU*1(**r***N*)

 0

 *U*2 1 0 − 1 3! *<sup>β</sup>*3*λ*<sup>3</sup>

<sup>−</sup>*βλU*1(**r***N*)

 *U*2 1 0 + 1 <sup>3</sup> �*U*1� 3 0 

*c*<sup>1</sup> = �*U*1�<sup>0</sup> , (33)

+ 2 �*U*1�

<sup>0</sup> <sup>=</sup> <sup>−</sup>*λβc*<sup>1</sup> <sup>+</sup> *<sup>λ</sup>*2*β*2*c*<sup>2</sup> <sup>−</sup> *<sup>λ</sup>*3*β*3*c*<sup>3</sup> <sup>+</sup> ...

*<sup>F</sup>* <sup>=</sup> *<sup>F</sup>*<sup>0</sup> <sup>+</sup> *<sup>F</sup>*<sup>1</sup> <sup>=</sup> *<sup>F</sup>*<sup>0</sup> <sup>+</sup> *<sup>λ</sup>c*<sup>1</sup> <sup>−</sup> *<sup>λ</sup>*2*βc*<sup>2</sup> <sup>+</sup> *<sup>λ</sup>*3*β*2*c*<sup>3</sup> <sup>+</sup> ..., (37)

∞ ∑ *n*=1

<sup>0</sup> <sup>=</sup> <sup>−</sup>*λβ*

3 0 

<sup>0</sup>, calculated on the reference system, is expressed as:

, (34)

 0

. (30)

+ ... (32)

. (31)

<sup>0</sup> in Taylor series. After

, etc. (35)

*cn*(−*λβ*)*n*−1, (36)

<sup>+</sup> *<sup>β</sup>*4*λ*<sup>4</sup> [...] <sup>−</sup> ...

 *U*3 1 0

(cf. footnote 1), this one reads:

statistical average

 exp 

ln exp 

Now if we set:

we find that:

*approximation*.

form:

− *βF* = ln

exp

<sup>−</sup>*βλU*1(**r***N*)

simplifications, equation (31) reduces to:

<sup>−</sup>*βλU*1(**r***N*)

*Z*(0) *<sup>N</sup>* (*V*, *T*) *<sup>N</sup>*!Λ3*<sup>N</sup>* <sup>+</sup> ln

and the second term represents the free energy of the perturbation:

<sup>−</sup>*βλU*1(**r***N*)

<sup>0</sup> <sup>=</sup> <sup>−</sup>*βλ* �*U*1�<sup>0</sup> <sup>+</sup>

 1 3! *U*3 1 0 − 1 <sup>2</sup> �*U*1�<sup>0</sup>

<sup>−</sup>*β*3*λ*<sup>3</sup>

<sup>−</sup>*βλU*1(**r***N*)

 exp 

possible to write equation (31) by expanding ln

*<sup>c</sup>*<sup>2</sup> <sup>=</sup> <sup>1</sup> 2! *U*2 1 0 − �*U*1� 2 0 

*<sup>c</sup>*<sup>3</sup> <sup>=</sup> <sup>1</sup> 3! *U*3 1 0

− *βF*<sup>1</sup> = ln

ln exp 

(36) into equation (30), as follows:

− *βF*<sup>1</sup> = ln

<sup>0</sup> <sup>=</sup> <sup>1</sup> <sup>−</sup> *βλ* �*U*1�<sup>0</sup> <sup>+</sup>

As a first application of thermodynamic perturbation method, we search the phenomenological van der Waals equation of state. In view of this, consider equation (37) at zero order in *β*. The simplest assumption to determine *c*<sup>1</sup> is to admit that the total potential energy may be decomposed into a sum of pair potentials in the form:

$$\mathcal{U}(\mathbf{r}^N) = \mathcal{U}\_0(\mathbf{r}^N) + \mathcal{U}\_1(\mathbf{r}^N) = \sum\_{i} \sum\_{j>i} \mu\_0(r\_{ij}) + \sum\_{i} \sum\_{j>i} \mu\_1(r\_{ij}).$$

Therefore, the free energy of the perturbation to zero order in *β* is given by equation (33), that is to say:

$$c\_1 = \langle \mathcal{U}\_1 \rangle\_0 = \frac{\int \int\_{\Im \mathcal{N}} \left[ \sum\_i \sum\_{j>i} u\_1(r\_{ij}) \right] \exp \left[ -\beta \mathcal{U}\_0(\mathbf{r}^N) \right] d\mathbf{r}^N}{\int \int\_{\Im \mathcal{N}} \exp \left[ -\beta \mathcal{U}\_0(\mathbf{r}^N) \right] d\mathbf{r}^N}.$$

To simplify the above relation, we proceed as for calculating the internal energy of liquids (Eq. 8) by revealing the pair correlation function of the reference system in the numerator. If we assume that the sum of pair potentials is composed of equivalent terms equal to <sup>∑</sup>*<sup>i</sup>* <sup>∑</sup>*j*>*<sup>i</sup> <sup>u</sup>*1(*rij*) = *<sup>N</sup>*(*N*−1) <sup>2</sup> *u*1(*r*12), the expression of *c*<sup>1</sup> is simplified as:

$$\mathcal{L}\_1 = \frac{1}{Z\_0(V, T)} \int \int \frac{\mathcal{N}(N - 1)}{2} \mu\_1(r\_{12}) \left\{ \int \int\_{\Im(N - 2)} \exp\left[ -\beta l l\_0(\mathbf{r}^N) \right] d\mathbf{r}\_3 \dots d\mathbf{r}\_N \right\} d\mathbf{r}\_1 d\mathbf{r}\_2 \dots$$

The integral in between the braces is then expressed as a function of the pair correlation function (cf. footnote 2), and *c*<sup>1</sup> reduces to:

$$\mathbf{c}\_{1} = \frac{\rho^{2}}{2} \int d\mathbf{R} \int u\_{1}(r)g\_{0}(r)d\mathbf{r}.\tag{38}$$

Yet, to find the equation of van der Waals we have to choose the hard-sphere system of diameter *σ*, as reference system, and suppose that the perturbation is a long-range potential, weakly attractive, the form of which is not useful to specify (Fig. 6a). Since one was unaware of the existence of the pair correlation function when the model was developed by van der Waals, it is reasonable to estimate *g*0(*r*) by a function equal to zero within the particle, and to one at the outside. According to van der Waals, suppose further that the available volume per particle8 is *b* = <sup>2</sup> <sup>3</sup>*πσ*<sup>3</sup> and the unoccupied volume is (*<sup>V</sup>* <sup>−</sup> *Nb*).

With these simplifications in mind, the configuration integral and free energy of the reference system are respectively (cf. footnote 1):

$$\begin{aligned} Z\_0(V, T) &= \int \int\_{3N} \exp\left[-\beta l l\_0(\mathbf{r}^N)\right] d\mathbf{r}^N = (V - Nb)^N, \\ \text{and} & \qquad F\_0 = -k\_B T \ln\left[\frac{(V - Nb)^N}{N! \Lambda^{3N}}\right] = -Nk\_B T \left[\ln\frac{(V - Nb)}{N} - 3\ln\Lambda + 1\right]. \end{aligned}$$

<sup>8</sup> The parameter *b* introduced by van der Waals is the covolume. Its expression comes from the fact that if two particles are in contact, half of the excluded volume <sup>4</sup> <sup>3</sup>*πσ*<sup>3</sup> must be assigned to each particle (Fig. 6b).

of Simple Liquids 23

Thermodynamic Perturbation Theory of Simple Liquids 861

*i*

before substituting it in the configuration integral. The advantage of this method is to calculate the coefficients *cn* and free energy (Eq. 37), not with the mean values of the perturbation, but with the fluctuation number of particles. Thus, each perturbating potential *u*1(*ri*) is constant

> *NiNj* 0

With the *local compressibility* approximation (LC), where *ρ* and *g*<sup>0</sup> depend on *p*, the expression of *c*<sup>2</sup> obtained by Barker and Henderson according to the method described above is written:

Incidentally, note the *macroscopic compressibility* approximation (MC), where only *ρ* is assumed to be dependent on *p*, has also been tested on a system made of the hard-sphere reference system and the square-well potential as perturbation. At low densities, the results of both approximations are comparable. But at intermediate densities, the results obtained with the LC approximation are in better agreement with the simulation results than the MC approximation. Note also that the coefficient *c*<sup>3</sup> has been calculated by Mansoori and Canfield

At this stage of the presentation of the thermodynamic perturbation theory, we are in position to calculate the first terms of the development of the free energy *F* (Eq. 37), using the hard-sphere system as reference system. But there is not yet a criterion for choosing the diameter *d* of hard spheres. This point is important because all potentials have a repulsive part that must be replaced by a hard-sphere potential of diameter properly chosen. Decisive progress has been made to solve this problem in three separate ways followed, respectively, by Barker and Henderson (19), Mansoori and Canfield (20) and Week, Chandler and Andersen

**Prescription of Barker and Henderson**. To choose the best reference system, that is to say, the optimal diameter of hard spheres, Barker and Henderson (19) proposed to replace the potential separation *u*(*r*) = *u*0(*r*) + *λu*1(*r*), where *u*0(*r*) is the reference potential, *u*1(*r*) the perturbation potential and *λ* the coupling parameter, by a more complicated separation

exp

where Ξ(*x*) is the Heaviside function, which is zero when *x* < 0 and is worth one when *x* > 0. Note that here *σ* is the value of *r* at which the real potential *u*(*r*) vanishes and *d* is the hard-sphere diameter of the reference potential, to be determined. Moreover, the parameters *λ* and *α* are coupling parameters that are 0 or 1. If one looks at equation (41) at the same time as figure(7a), it is seen that the function *v*(*r*) reduces to the real potential *u*(*r*) when *α* = *λ* = 1,

−*βu*(*d* +

*r* − *d <sup>α</sup>* ) 

+ Ξ (*r* − *σ*) {exp [−*βλu*(*r*)] − 1} , (41)

*r* − *d <sup>α</sup>* <sup>−</sup> *<sup>σ</sup>*

*r* − *d <sup>α</sup>* <sup>−</sup> *<sup>σ</sup>* *Niu*1(*ri*),

2

− �*Ni*�<sup>0</sup>

<sup>1</sup>(*r*)*g*0(*r*)*r*

*ρu*<sup>2</sup>

 *Nj* 0 

<sup>2</sup>*dr*

<sup>0</sup> = ∑*<sup>i</sup>* ∑*j*�*Ni*�0�*Nj*�0*u*1(*ri*)*u*1(*rj*). In

*u*1(*ri*)*u*1(*rj*).

. (40)

*<sup>U</sup>*1(**r***N*) = ∑

in the interval which it belongs, so that we can write �*U*1�

view of this, the coefficient *c*<sup>2</sup> defined by equation (34) is:

*<sup>c</sup>*2(*LC*) = *<sup>π</sup>N<sup>ρ</sup>*

(18) with the macroscopic compressibility approximation.

associated with a potential *v*(*r*) whose the Boltzmann factor is:

+ Ξ *d* +

 1 − Ξ *d* +

exp [−*βv*(*r*)] =

*β*

 *∂ρ ∂p* 0 *∂ ∂ρ*

*<sup>c</sup>*<sup>2</sup> <sup>=</sup> <sup>1</sup> 2!

(21).

*<sup>U</sup>*<sup>2</sup> 1 0 − �*U*1� 2 0 <sup>=</sup> <sup>1</sup> 2! ∑ *i* ∑ *j*

Fig. 6. Schematic representation of the pair potential by a hard-sphere potentiel plus a perturbation. (b) Definition of the covolume by the quantity *b* = <sup>1</sup> 2 <sup>4</sup> 3*πσ*<sup>3</sup> .

As for the coefficient *c*<sup>1</sup> (Eq. 38), it is simplified as:

$$\begin{split} \mathfrak{c}\_{1} &= 2\pi \rho^{2} V \int\_{\sigma}^{\infty} u\_{1}(r) r^{2} dr = -a\rho N, \\ \text{with} \qquad a &= -2\pi \int\_{\sigma}^{\infty} u\_{1}(r) r^{2} dr. \end{split} \tag{39}$$

Therefore, the expression of free energy (Eq. 37) corresponding to the model of van der Waals is:

$$F = F\_0 + \mathcal{c}\_1 = -Nk\_B T \left[ \ln \frac{(V - Nb)}{N} - 3 \ln \Lambda + 1 \right] - a\rho N\_\prime \bar{\lambda}$$

and the van der Waals equation of state reduces to:

$$p = -\left(\frac{\partial F}{\partial V}\right)\_T = \frac{Nk\_B T}{V - Nb} - a\frac{N}{V^2}^2.$$

With *b* = *a* = 0 in the previous equation, it is obvious that one recovers the equation of state of ideal gas. In return, if one wishes to improve the quality of the van der Waals equation of state, one may use the expression of the free energy (Eq. 25) and pair correlation function *gHS*(*r*) of the hard-sphere system to calculate the value of the parameter *a* with equation (38). Another way to improve performance is to calculate the term *c*2. Precisely what will be done in the next subsection.

#### **5.3 Method of Barker and Henderson**

To evaluate the mean values of the perturbation *U*1(**r***N*) in equations (34) and (35), Barker and Henderson (17) suggested to discretize the domain of interatomic distances into sufficiently small intervals (*r*1,*r*2),(*r*2,*r*3), ...,(*ri*,*ri*<sup>+</sup>1),..., and assimilate the perturbating elemental potential in each interval by a constant. Assuming that the perturbating potential in the interval (*ri*,*ri*<sup>+</sup>1) is *u*1(*ri*) and that the number of atoms subjected to this potential is *Ni*, the total perturbation can be written as the sum of elemental potentials:

22 Thermodynamics book 1

Fig. 6. Schematic representation of the pair potential by a hard-sphere potentiel plus a

 ∞ *σ*

Therefore, the expression of free energy (Eq. 37) corresponding to the model of van der Waals

ln (*<sup>V</sup>* <sup>−</sup> *Nb*)

<sup>=</sup> *NkBT <sup>V</sup>* <sup>−</sup> *Nb* <sup>−</sup> *<sup>a</sup>*

With *b* = *a* = 0 in the previous equation, it is obvious that one recovers the equation of state of ideal gas. In return, if one wishes to improve the quality of the van der Waals equation of state, one may use the expression of the free energy (Eq. 25) and pair correlation function *gHS*(*r*) of the hard-sphere system to calculate the value of the parameter *a* with equation (38). Another way to improve performance is to calculate the term *c*2. Precisely what will be done

To evaluate the mean values of the perturbation *U*1(**r***N*) in equations (34) and (35), Barker and Henderson (17) suggested to discretize the domain of interatomic distances into sufficiently small intervals (*r*1,*r*2),(*r*2,*r*3), ...,(*ri*,*ri*<sup>+</sup>1),..., and assimilate the perturbating elemental potential in each interval by a constant. Assuming that the perturbating potential in the interval (*ri*,*ri*<sup>+</sup>1) is *u*1(*ri*) and that the number of atoms subjected to this potential is *Ni*, the

 ∞ *σ*

 *∂F ∂V T* *u*1(*r*)*r*

<sup>2</sup>*dr*.

*<sup>N</sup>* <sup>−</sup> 3 ln <sup>Λ</sup> <sup>+</sup> <sup>1</sup>

*N V*<sup>2</sup> 2 .

*u*1(*r*)*r*

*c*<sup>1</sup> = 2*πρ*2*V*

2 <sup>4</sup> 3*πσ*<sup>3</sup> .

− *aρN*,

<sup>2</sup>*dr* <sup>=</sup> <sup>−</sup>*aρN*, (39)

perturbation. (b) Definition of the covolume by the quantity *b* = <sup>1</sup>

with *a* = −2*π*

*F* = *F*<sup>0</sup> + *c*<sup>1</sup> = −*NkBT*

*p* = −

total perturbation can be written as the sum of elemental potentials:

and the van der Waals equation of state reduces to:

As for the coefficient *c*<sup>1</sup> (Eq. 38), it is simplified as:

is:

in the next subsection.

**5.3 Method of Barker and Henderson**

$$\mathcal{U}\_1(\mathbf{r}^N) = \sum\_i \mathcal{N}\_i \mu\_1(r\_i)\_{\prime}$$

before substituting it in the configuration integral. The advantage of this method is to calculate the coefficients *cn* and free energy (Eq. 37), not with the mean values of the perturbation, but with the fluctuation number of particles. Thus, each perturbating potential *u*1(*ri*) is constant in the interval which it belongs, so that we can write �*U*1� 2 <sup>0</sup> = ∑*<sup>i</sup>* ∑*j*�*Ni*�0�*Nj*�0*u*1(*ri*)*u*1(*rj*). In view of this, the coefficient *c*<sup>2</sup> defined by equation (34) is:

$$\mathcal{L}\_2 = \frac{1}{2!} \left[ \left< \mathcal{U}\_1^2 \right>\_0 - \left< \mathcal{U}\_1 \right>\_0^2 \right] = \frac{1}{2!} \sum\_{i} \sum\_{j} \left[ \left< \mathcal{N}\_i \mathcal{N}\_j \right>\_0 - \left< \mathcal{N}\_i \right>\_0 \left< \mathcal{N}\_j \right>\_0 \right] u\_1(r\_i) u\_1(r\_j).$$

With the *local compressibility* approximation (LC), where *ρ* and *g*<sup>0</sup> depend on *p*, the expression of *c*<sup>2</sup> obtained by Barker and Henderson according to the method described above is written:

$$c\_2(L\mathbb{C}) = \frac{\pi N \rho}{\beta} \left(\frac{\partial \rho}{\partial p}\right)\_0 \frac{\partial}{\partial \rho} \left[\int \rho u\_1^2(r) g\_0(r) r^2 dr\right]. \tag{40}$$

Incidentally, note the *macroscopic compressibility* approximation (MC), where only *ρ* is assumed to be dependent on *p*, has also been tested on a system made of the hard-sphere reference system and the square-well potential as perturbation. At low densities, the results of both approximations are comparable. But at intermediate densities, the results obtained with the LC approximation are in better agreement with the simulation results than the MC approximation. Note also that the coefficient *c*<sup>3</sup> has been calculated by Mansoori and Canfield (18) with the macroscopic compressibility approximation.

At this stage of the presentation of the thermodynamic perturbation theory, we are in position to calculate the first terms of the development of the free energy *F* (Eq. 37), using the hard-sphere system as reference system. But there is not yet a criterion for choosing the diameter *d* of hard spheres. This point is important because all potentials have a repulsive part that must be replaced by a hard-sphere potential of diameter properly chosen. Decisive progress has been made to solve this problem in three separate ways followed, respectively, by Barker and Henderson (19), Mansoori and Canfield (20) and Week, Chandler and Andersen (21).

**Prescription of Barker and Henderson**. To choose the best reference system, that is to say, the optimal diameter of hard spheres, Barker and Henderson (19) proposed to replace the potential separation *u*(*r*) = *u*0(*r*) + *λu*1(*r*), where *u*0(*r*) is the reference potential, *u*1(*r*) the perturbation potential and *λ* the coupling parameter, by a more complicated separation associated with a potential *v*(*r*) whose the Boltzmann factor is:

$$\begin{split} \mathbb{E}\exp\left[-\beta v(r)\right] &= \left[1 - \Xi\left(d + \frac{r-d}{a} - \sigma\right)\right] \exp\left[-\beta u(d + \frac{r-d}{a})\right] \\ &+ \Xi\left(d + \frac{r-d}{a} - \sigma\right) + \Xi\left(r - \sigma\right) \left\{\exp\left[-\beta \lambda u(r)\right] - 1\right\}, \end{split} \tag{41}$$

where Ξ(*x*) is the Heaviside function, which is zero when *x* < 0 and is worth one when *x* > 0. Note that here *σ* is the value of *r* at which the real potential *u*(*r*) vanishes and *d* is the hard-sphere diameter of the reference potential, to be determined. Moreover, the parameters *λ* and *α* are coupling parameters that are 0 or 1. If one looks at equation (41) at the same time as figure(7a), it is seen that the function *v*(*r*) reduces to the real potential *u*(*r*) when *α* = *λ* = 1,

of Simple Liquids 25

Thermodynamic Perturbation Theory of Simple Liquids 863

Therefore, using equation (38) to evaluate *c*<sup>1</sup> and equation (40) to evaluate *c*2, the expression

*u*1(*r*)*gHS*(*η*;*r*)*r*

*∂ ∂ρ* <sup>∞</sup> *d ρu*<sup>2</sup>

where the first term on the RHS represents the free energy of the hard-sphere system (Eq. 25),

state (Eq. 24). Recall that the pair correlation function of hard-sphere system, *gHS*(*η*;*r*), must be only determined numerically. It depends on the density *ρ* and diameter *d* via the packing

An important consequence of the high temperature approximation to first order in *β* is to mark out the free energy of the real system by an upper limit that can not be exceeded, because the sum of the terms beyond *c*<sup>1</sup> is always negative. The easiest way to proof this, is to consider the expression of the free energy (Eq. 30) and to write the perturbation *U*1(**r***N*) around its mean

*U*<sup>1</sup> = �*U*1�<sup>0</sup> + Δ*U*1.

However, considering the series expansion of an exponential, the above relation is

A thorough study of this inequality shows that it is always valid, and it is unnecessary to consider values of *n* greater than zero. Equation (46), at the base of the variational method, allows us to find the value of the parameter *d* that makes the free energy *F minimum*. If the reference system is that of hard spheres, the value of the free energy obtained with this value

<sup>6</sup> *<sup>ρ</sup>d*3) is considered as the best estimate of the free energy of the real system. Its

*<sup>k</sup>*! (*<sup>n</sup>* <sup>=</sup> 0, 1, 2, ...).

= ln �1 − *β*Δ*U*1�<sup>0</sup> � −*β* �Δ*U*1�<sup>0</sup> = 0,

<sup>−</sup> *<sup>β</sup><sup>F</sup>* <sup>=</sup> <sup>−</sup>*βF*<sup>0</sup> <sup>−</sup> *<sup>β</sup>* �*U*1�<sup>0</sup> <sup>+</sup> ln

*<sup>F</sup>* <sup>≤</sup> *FHS* <sup>+</sup> *<sup>ρ</sup><sup>N</sup>*

<sup>9</sup> The Mac-Laurin series of the exponential naturally leads to the inequality:

exp (−*β*Δ*U*1) ≥

(−*β*Δ*U*1) *k*

*k*!

For *n* = 0, the last term of equation (45) behaves as:

ln <sup>1</sup> ∑ *k*=0

since the mean value of the deviation, �Δ*U*1�<sup>0</sup> , is zero.

2 

2*n*+1 ∑ *k*=0

0

(−*β*Δ*U*1) *k*

<sup>6</sup> *<sup>ρ</sup>d*3). Since *gHS*(*η*;*r*) = 0 when *<sup>r</sup>* <sup>&</sup>lt; *<sup>d</sup>*, either 0 or *<sup>d</sup>* can be used as lower limit

<sup>2</sup>*dr*

<sup>1</sup>(*r*)*gHS*(*η*;*r*)*r*

*HS* can be deduced from the Carnahan and Starling equation of

exp [−*β*Δ*U*1]

0 

*F* ≤ *F*<sup>0</sup> + �*U*1�<sup>0</sup> . (46)

*u*1(*r*)*gHS*(*η*;*r*)*d***r**, (47)

<sup>2</sup>*dr*

, (44)

. (45)

 ∞ *d*

 *∂ρ ∂p* 

*HS*

of the free energy *F* of the real system (Eq. 42) is:

and the partial derivative *∂ρ*

of integration in equation (44).

**5.4 Prescription of Mansoori and Canfield.**

After replacing *U*<sup>1</sup> in equation (30), we obtain:

transformed into the so-called *Gibbs-Bogoliubov inequality*9:

fraction *η*(= *<sup>π</sup>*

value �*U*1�<sup>0</sup> as:

of *d* (or *η* = *<sup>π</sup>*

expression is:

*F* = *FHS* + 2*πρN*

− *πNρ*

*∂p* 

Fig. 7. Separation of the potential *u*(*r*) according to (a) the method of Barker and Henderson and (b) the method of Weeks, Chandler and Andersen.

and it behaves approximately as the hard-sphere potential of diameter *d* when *α* ∼ *λ* ∼ 0. The substitution of equation (41) in the configuration integral (Eq. 29), followed by the related calculations not reproduced here, enable us to express the free energy *F* of the real system as a series expansion in powers of *α* and *λ*, which makes the generalization of equation (28), namely:

$$F = F\_{HS} + \lambda \left(\frac{\partial F}{\partial \lambda}\right) + a \left(\frac{\partial F}{\partial a}\right) + \frac{\lambda^2}{2} \left(\frac{\partial^2 F}{\partial \lambda^2}\right) + \frac{a^2}{2} \left(\frac{\partial^2 F}{\partial a^2}\right) + \dots \tag{42}$$

By comparing equations (37) and (42), we see that the first derivative *<sup>∂</sup><sup>F</sup> ∂λ* coincides with *c*<sup>1</sup> and the second derivative <sup>1</sup> 2 *<sup>∂</sup>*2*<sup>F</sup> ∂λ*<sup>2</sup> with (−*βc*2). Concerning the derivatives of *F* with respect to *α*, they are complicated functions of the pair potential and the pair correlation function of the hard-sphere system. The first derivative *<sup>∂</sup><sup>F</sup> ∂α* , whose the explicit form given without proof, reads:

$$\left(\frac{\partial F}{\partial \boldsymbol{a}}\right) = -2\pi N \rho k\_B T d^2 g^{HS}(d) \left[d - \int\_0^{\sigma} \left\{1 - \exp\left[-\beta \mu(\boldsymbol{r})\right]\right\} d\boldsymbol{r}\right].$$

Since the Barker and Henderson prescription is based on the proposal to cancel the term *<sup>∂</sup><sup>F</sup> ∂α* , the criterion for choosing the hard-sphere diameter *d* is reduced to the following equation:

$$d = \int\_0^{\sigma} \left\{ 1 - \exp \left[ -\beta u(r) \right] \right\} dr. \tag{43}$$

In applying this criterion to the Lennard-Jones potential, it is seen that *d* depends on temperature but not on the density. Also, the calculations show that the terms of the expansion of *F* in *α*<sup>2</sup> and *αλ* are negligible compared to the term in *λ*2.

24 Thermodynamics book 1

Fig. 7. Separation of the potential *u*(*r*) according to (a) the method of Barker and Henderson

and it behaves approximately as the hard-sphere potential of diameter *d* when *α* ∼ *λ* ∼ 0. The substitution of equation (41) in the configuration integral (Eq. 29), followed by the related calculations not reproduced here, enable us to express the free energy *F* of the real system as a series expansion in powers of *α* and *λ*, which makes the generalization of equation (28),

to *α*, they are complicated functions of the pair potential and the pair correlation function

 *d* − *σ*

Since the Barker and Henderson prescription is based on the proposal to cancel the term

the criterion for choosing the hard-sphere diameter *d* is reduced to the following equation:

In applying this criterion to the Lennard-Jones potential, it is seen that *d* depends on temperature but not on the density. Also, the calculations show that the terms of the expansion

 *<sup>∂</sup><sup>F</sup> ∂α* 

 *∂*2*F ∂λ*<sup>2</sup> + *α*2 2

 *∂*2*F ∂α*<sup>2</sup> 

with (−*βc*2). Concerning the derivatives of *F* with respect

<sup>0</sup> {<sup>1</sup> <sup>−</sup> exp [−*βu*(*r*)]} *dr*

<sup>0</sup> {<sup>1</sup> <sup>−</sup> exp [−*βu*(*r*)]} *dr*. (43)

 *<sup>∂</sup><sup>F</sup> ∂λ* 

, whose the explicit form given without

 .

+ ... (42)

coincides with *c*<sup>1</sup>

 *<sup>∂</sup><sup>F</sup> ∂α* ,

 *∂F ∂α* + *λ*2 2

and (b) the method of Weeks, Chandler and Andersen.

 *∂F ∂λ* + *α*

2 *<sup>∂</sup>*2*<sup>F</sup> ∂λ*<sup>2</sup> 

of the hard-sphere system. The first derivative

By comparing equations (37) and (42), we see that the first derivative

<sup>=</sup> <sup>−</sup>2*πNρkBTd*<sup>2</sup>*gHS*(*d*)

*d* = *σ*

of *F* in *α*<sup>2</sup> and *αλ* are negligible compared to the term in *λ*2.

*F* = *FHS* + *λ*

and the second derivative <sup>1</sup>

 *∂F ∂α* 

namely:

proof, reads:

Therefore, using equation (38) to evaluate *c*<sup>1</sup> and equation (40) to evaluate *c*2, the expression of the free energy *F* of the real system (Eq. 42) is:

$$\begin{split} F &= F\_{\rm HS} + 2\pi\rho \text{N} \int\_{d}^{\infty} u\_{1}(r) g\_{\rm HS}(\eta; r) r^{2} dr \\ &\quad - \pi \text{N} \rho \left( \frac{\partial \rho}{\partial p} \right)\_{\rm HS} \frac{\partial}{\partial \rho} \left[ \int\_{d}^{\infty} \rho u\_{1}^{2}(r) g\_{\rm HS}(\eta; r) r^{2} dr \right], \end{split} \tag{44}$$

where the first term on the RHS represents the free energy of the hard-sphere system (Eq. 25), and the partial derivative *∂ρ ∂p HS* can be deduced from the Carnahan and Starling equation of state (Eq. 24). Recall that the pair correlation function of hard-sphere system, *gHS*(*η*;*r*), must be only determined numerically. It depends on the density *ρ* and diameter *d* via the packing fraction *η*(= *<sup>π</sup>* <sup>6</sup> *<sup>ρ</sup>d*3). Since *gHS*(*η*;*r*) = 0 when *<sup>r</sup>* <sup>&</sup>lt; *<sup>d</sup>*, either 0 or *<sup>d</sup>* can be used as lower limit of integration in equation (44).

#### **5.4 Prescription of Mansoori and Canfield.**

An important consequence of the high temperature approximation to first order in *β* is to mark out the free energy of the real system by an upper limit that can not be exceeded, because the sum of the terms beyond *c*<sup>1</sup> is always negative. The easiest way to proof this, is to consider the expression of the free energy (Eq. 30) and to write the perturbation *U*1(**r***N*) around its mean value �*U*1�<sup>0</sup> as:

$$\mathcal{U}\_1 = \langle \mathcal{U}\_1 \rangle\_0 + \Delta \mathcal{U}\_1$$

After replacing *U*<sup>1</sup> in equation (30), we obtain:

$$-\beta F = -\beta F\_0 - \beta \left< \mathcal{U}\_1 \right>\_0 + \ln \left< \exp \left[ -\beta \Delta \mathcal{U}\_1 \right]\_0 \right>. \tag{45}$$

However, considering the series expansion of an exponential, the above relation is transformed into the so-called *Gibbs-Bogoliubov inequality*9:

$$F \le F\_0 + \langle \mathcal{U}\_1 \rangle\_0. \tag{46}$$

A thorough study of this inequality shows that it is always valid, and it is unnecessary to consider values of *n* greater than zero. Equation (46), at the base of the variational method, allows us to find the value of the parameter *d* that makes the free energy *F minimum*. If the reference system is that of hard spheres, the value of the free energy obtained with this value of *d* (or *η* = *<sup>π</sup>* <sup>6</sup> *<sup>ρ</sup>d*3) is considered as the best estimate of the free energy of the real system. Its expression is:

$$F \le F\_{\rm HS} + \frac{\rho \mathcal{N}}{2} \int u\_1(r) g\_{\rm HS}(\eta; r) d\mathbf{r}\_\prime \tag{47}$$

$$\exp\left(-\beta\Delta U I\_1\right) \ge \sum\_{k=0}^{2n+1} \frac{\left(-\beta\Delta U I\_1\right)^k}{k!} \qquad \text{(\$n=0\$, 1, 2, ...\$)}\text{.}$$

For *n* = 0, the last term of equation (45) behaves as:

$$\ln \left\langle \sum\_{k=0}^{1} \frac{\left(-\beta \Delta U I\_1\right)^k}{k!} \right\rangle\_0 = \ln \left\langle 1 - \beta \Delta U I\_1 \right\rangle\_0 \cong -\beta \left\langle \Delta U I\_1 \right\rangle\_0 = 0,$$

since the mean value of the deviation, �Δ*U*1�<sup>0</sup> , is zero.

<sup>9</sup> The Mac-Laurin series of the exponential naturally leads to the inequality:

of Simple Liquids 27

Thermodynamic Perturbation Theory of Simple Liquids 865

It should be stressed that the cavity function depends only weakly on the potential, this is why Weeks, Chander and Anderson (WCA) suggested to use the same cavity function *yHS*(*η*;*r*) for all potentials, and to express the pair correlation functions of each potential as a function of *yHS*(*η*;*r*). As a result, the pair correlation function *g*0(*r*) related to the repulsive contribution

The suggestion of WCA for choosing the hard-sphere diameter *<sup>d</sup>* is to cancel *<sup>∂</sup><sup>F</sup>*

<sup>0</sup> {exp [−*βu*0(*r*)] <sup>−</sup> exp [−*βuHS*(*d*)]} *yHS*(*η*;*r*)*<sup>r</sup>*

<sup>2</sup>*dr* =

It is interesting to note that equation (53) has a precise physical meaning that appears by writing it in terms of the pair correlation functions *gHS*(*η*;*r*) and *g*0(*r*) drawn, respectively,

> ∞ 0

According to equation (7), the previous expression is equivalent to the equality *S*0(0) = *SHS*(0), meaning the equality between the isothermal compressibility of the repulsive

Ultimately, the expression of free energy (Eq. 28) of the real system is the sum of *FHS*,

with equation (38) in using the pair correlation function *g*0(*r*) � exp [−*βu*0(*r*)] *yHS*(*η*;*r*), i. e.:

Note that the value of *d* obtained with the WCA prescription (*blip function* method) is significantly larger than that obtained with the Barker and Henderson (BH) prescription. By the fact that *yHS*(*η*;*r*) depends on density, the value of *d* calculated with equation (53) depends on temperature and density, while that calculated with equation (43) depends only

In order to compare the relative merits of equations (44) and (54), simply note that the WCA prescription shows that not only the terms of the expansion of *F* in *α*<sup>2</sup> and in *αλ* are negligible, but also the term in *λ*2. This means that the WCA treatment is a theory of first order in *λ*, while the BH one is a theory of second order in *λ*. In addition, the BH treatment is coherent as *FHS*, *c*<sup>1</sup> and *c*<sup>2</sup> are calculated with the same reference system, whereas the WCA treatment is not coherent because it uses the hard-sphere system to calculate *FHS* and the repulsive potential *u*0(*r*) to calculate *c*1, via the pair correlation function *g*0(*r*). In contrast, the WCA treatment does not require the calculation of *c*<sup>2</sup> and improves the convergence of calculations. As an additional advantage, the WCA treatment can be used to predict the structure factor of simple

[*gHS*(*η*;*r*) − 1]*r*

*u*1(*r*) {exp [−*βu*0(*r*)] *yHS*(*η*;*r*)} *r*

<sup>2</sup>*dr*.

*∂λ*

= *c*1, calculated

<sup>2</sup>*dr*. (54)

amounts to solving the nonlinear equation given by equation (50), i. e.:

[*g*0(*r*) − 1]*r*

calculated with the value of *<sup>d</sup>* issued from equation (53), and the term *<sup>∂</sup><sup>F</sup>*

 ∞ 0

potential *u*0(*r*) and that of the hard-sphere potential *uHS*(*d*).

*g*0(*r*) � exp [−*βu*0(*r*)] *yHS*(*η*;*r*). (52)

*∂α*

<sup>2</sup>*dr* = 0. (53)

, which

*u*0(*r*) of the potential can be approximately written as:

∞

from equations (51) and (52). It reads:

on temperature.

liquids (22).

 ∞ 0

*F* = *FHS* + 2*πρN*

where *F*<sup>0</sup> = *FHS* is given by equation (25) and �*U*1�<sup>0</sup> = *c*<sup>1</sup> by equation (38). Since *gHS*(*η*;*r*) = 0 when *r* < *d* and *u*1(*r*) = *u*(*r*) − *uHS*(*r*) = *u*(*r*) when *r* ≥ *d*, equation (47) can be written as:

$$F \le F\_{\rm HS} + 2\pi\rho N \int\_0^\infty u(r) g\_{\rm HS}(\eta; r) r^2 dr. \tag{48}$$

Practically, we vary the value of *d* (or *η*) until the result of the integral is minimum. And the value of *F* thus obtained is the best estimate of the free energy of the real system, in the sens of the variational method (20).

#### **5.5 Prescription of Weeks, Chandler and Andersen.**

At the same time that the perturbation theory was developing, Weeks, Chandler and Andersen (21) formulated another prescription for finding the hard-sphere diameter *d*. Its originality lies in the idea that a particle in the liquid is less sensitive to the sign of the potential *u*(*r*) than to the sign of the strength, that is to say, to the derivative of the potential <sup>−</sup>*∂<sup>u</sup> ∂r ρ*,*T* . That is why the authors proposed to separate the real potential into a purely repulsive contribution and a purely attractive perturbation. This separation, shown in figure (7b), is defined by the relation *u*(*r*) = *u*0(*r*) + *u*1(*r*) with:

$$\begin{aligned} u\_0(r) &= \begin{cases} u(r) + \varepsilon & \text{if} \quad r < r\_{\text{min}} \\ 0 & \text{if} \quad r \ge r\_{\text{min}} \end{cases} \\ u\_1(r) &= \begin{cases} -\varepsilon & \text{if} \quad r < r\_{\text{min}} \\ u(r) & \text{if} \quad r \ge r\_{\text{min}} \end{cases} \end{aligned}$$

where *ε* is the depth of the potential well, i. e. the value of *u*(*r*min) = −*ε*. To follow the same sketch that for the Barker and Henderson prescription, define the potential *v*(*r*) by the Boltzmann factor:

$$\exp\left[-\beta v(r)\right] = \exp\left[-\beta u\_{HS}(d)\right] + \mathfrak{a}\left\{\exp\left[-\beta\left(u\_0(r) + \lambda u\_1(r)\right)\right] - \exp\left[-\beta u\_{HS}(d)\right]\right\}.\tag{49}$$

However, from the simultaneous observation of equation (49) and figure (7b) one remarks that the potential *v*(*r*) reduces to the real potential *u*(*r*) when *α* = *λ* = 1, and behaves like the hard-sphere potential of diameter *d* when *α* = 0.

The substitution of equation (49) in the configuration integral, and subsequent calculations, enables us to express the free energy *F* of the system as a series expansion in powers of *α* and *λ* identical to that of equation (42). The first term of this expansion, *FHS*, is the free energy of the hard-sphere system. As for the first derivative of *F* with respect to *α*, it is written in this case:

$$\left(\frac{\partial F}{\partial \mathbf{a}}\right) = -\frac{2\pi N\rho}{\beta} \int\_0^\infty \left\{ \exp\left[-\beta u\_0(r)\right] - \exp\left[-\beta u\_{HS}(d)\right] \right\} y\_{HS}(\eta; r) r^2 dr,\tag{50}$$

where *yHS*(*η*;*r*) is the *cavity function* of great importance in the microscopic theory of liquids10. It is defined by means of the pair correlation function *gHS*(*η*;*r*) and the Boltzmann factor exp [−*βuHS*(*d*)] of the hard-sphere potential as:

$$\|\!\|y\_{HS}(\eta;r) = \exp\left[\beta\mu\_{HS}(d)\right]\|\_{HS}(\eta;r). \tag{51}$$

<sup>10</sup> The cavity function *yHS*(*η*;*r*) does not exist in analytical form. It must be calculated at all reduced densities (*ρd*3) by solving an integro-differential equation.

26 Thermodynamics book 1

where *F*<sup>0</sup> = *FHS* is given by equation (25) and �*U*1�<sup>0</sup> = *c*<sup>1</sup> by equation (38). Since *gHS*(*η*;*r*) = 0 when *r* < *d* and *u*1(*r*) = *u*(*r*) − *uHS*(*r*) = *u*(*r*) when *r* ≥ *d*, equation (47) can be written as:

> ∞ 0

Practically, we vary the value of *d* (or *η*) until the result of the integral is minimum. And the value of *F* thus obtained is the best estimate of the free energy of the real system, in the sens

At the same time that the perturbation theory was developing, Weeks, Chandler and Andersen (21) formulated another prescription for finding the hard-sphere diameter *d*. Its originality lies in the idea that a particle in the liquid is less sensitive to the sign of the potential *u*(*r*) than to the sign of the strength, that is to say, to the derivative of the potential

repulsive contribution and a purely attractive perturbation. This separation, shown in figure

. That is why the authors proposed to separate the real potential into a purely

 *u*(*r*) + *ε* if *r* < *r*min 0 if *r* ≥ *r*min

 <sup>−</sup>*<sup>ε</sup>* if *<sup>r</sup>* <sup>&</sup>lt; *<sup>r</sup>*min *u*(*r*) if *r* ≥ *r*min,

To follow the same sketch that for the Barker and Henderson prescription, define the potential

exp [−*βv*(*r*)] = exp [−*βuHS*(*d*)] + *α* {exp [−*β* (*u*0(*r*) + *λu*1(*r*))] − exp [−*βuHS*(*d*)]} . (49)

However, from the simultaneous observation of equation (49) and figure (7b) one remarks that the potential *v*(*r*) reduces to the real potential *u*(*r*) when *α* = *λ* = 1, and behaves like the

The substitution of equation (49) in the configuration integral, and subsequent calculations, enables us to express the free energy *F* of the system as a series expansion in powers of *α* and *λ* identical to that of equation (42). The first term of this expansion, *FHS*, is the free energy of the hard-sphere system. As for the first derivative of *F* with respect to *α*, it is written in this

where *yHS*(*η*;*r*) is the *cavity function* of great importance in the microscopic theory of liquids10. It is defined by means of the pair correlation function *gHS*(*η*;*r*) and the Boltzmann factor

<sup>10</sup> The cavity function *yHS*(*η*;*r*) does not exist in analytical form. It must be calculated at all reduced

<sup>0</sup> {exp [−*βu*0(*r*)] <sup>−</sup> exp [−*βuHS*(*d*)]} *yHS*(*η*;*r*)*<sup>r</sup>*

*yHS*(*η*;*r*) = exp [*βuHS*(*d*)] *gHS*(*η*;*r*). (51)

*u*(*r*)*gHS*(*η*;*r*)*r*

<sup>2</sup>*dr*. (48)

<sup>2</sup>*dr*, (50)

*F* ≤ *FHS* + 2*πρN*

of the variational method (20).

*v*(*r*) by the Boltzmann factor:

 *∂F ∂α* 

 <sup>−</sup>*∂<sup>u</sup> ∂r ρ*,*T*

case:

**5.5 Prescription of Weeks, Chandler and Andersen.**

(7b), is defined by the relation *u*(*r*) = *u*0(*r*) + *u*1(*r*) with:

hard-sphere potential of diameter *d* when *α* = 0.

<sup>=</sup> <sup>−</sup>2*πN<sup>ρ</sup> β*

exp [−*βuHS*(*d*)] of the hard-sphere potential as:

densities (*ρd*3) by solving an integro-differential equation.

∞

*u*0(*r*) =

*u*1(*r*) =

where *ε* is the depth of the potential well, i. e. the value of *u*(*r*min) = −*ε*.

It should be stressed that the cavity function depends only weakly on the potential, this is why Weeks, Chander and Anderson (WCA) suggested to use the same cavity function *yHS*(*η*;*r*) for all potentials, and to express the pair correlation functions of each potential as a function of *yHS*(*η*;*r*). As a result, the pair correlation function *g*0(*r*) related to the repulsive contribution *u*0(*r*) of the potential can be approximately written as:

$$\log\_0(r) \simeq \exp\left[-\beta u\_0(r)\right] y\_{HS}(\eta; r). \tag{52}$$

The suggestion of WCA for choosing the hard-sphere diameter *<sup>d</sup>* is to cancel *<sup>∂</sup><sup>F</sup> ∂α* , which amounts to solving the nonlinear equation given by equation (50), i. e.:

$$\int\_0^\infty \left\{ \exp\left[-\beta u\_0(r)\right] - \exp\left[-\beta u\_{HS}(d)\right] \right\} y\_{HS}(\eta; r) r^2 dr = 0. \tag{53}$$

It is interesting to note that equation (53) has a precise physical meaning that appears by writing it in terms of the pair correlation functions *gHS*(*η*;*r*) and *g*0(*r*) drawn, respectively, from equations (51) and (52). It reads:

$$\int\_0^\infty \left[g\_0(r) - 1\right] r^2 dr = \int\_0^\infty \left[g\_{HS}(\eta; r) - 1\right] r^2 dr.$$

According to equation (7), the previous expression is equivalent to the equality *S*0(0) = *SHS*(0), meaning the equality between the isothermal compressibility of the repulsive potential *u*0(*r*) and that of the hard-sphere potential *uHS*(*d*).

Ultimately, the expression of free energy (Eq. 28) of the real system is the sum of *FHS*, calculated with the value of *<sup>d</sup>* issued from equation (53), and the term *<sup>∂</sup><sup>F</sup> ∂λ* = *c*1, calculated with equation (38) in using the pair correlation function *g*0(*r*) � exp [−*βu*0(*r*)] *yHS*(*η*;*r*), i. e.:

$$F = F\_{HS} + 2\pi\rho N \int\_0^\infty u\_1(r) \left\{ \exp\left[-\beta u\_0(r)\right] y\_{HS}(\eta; r) \right\} r^2 dr. \tag{54}$$

Note that the value of *d* obtained with the WCA prescription (*blip function* method) is significantly larger than that obtained with the Barker and Henderson (BH) prescription. By the fact that *yHS*(*η*;*r*) depends on density, the value of *d* calculated with equation (53) depends on temperature and density, while that calculated with equation (43) depends only on temperature.

In order to compare the relative merits of equations (44) and (54), simply note that the WCA prescription shows that not only the terms of the expansion of *F* in *α*<sup>2</sup> and in *αλ* are negligible, but also the term in *λ*2. This means that the WCA treatment is a theory of first order in *λ*, while the BH one is a theory of second order in *λ*. In addition, the BH treatment is coherent as *FHS*, *c*<sup>1</sup> and *c*<sup>2</sup> are calculated with the same reference system, whereas the WCA treatment is not coherent because it uses the hard-sphere system to calculate *FHS* and the repulsive potential *u*0(*r*) to calculate *c*1, via the pair correlation function *g*0(*r*). In contrast, the WCA treatment does not require the calculation of *c*<sup>2</sup> and improves the convergence of calculations. As an additional advantage, the WCA treatment can be used to predict the structure factor of simple liquids (22).

of Simple Liquids 29

Thermodynamic Perturbation Theory of Simple Liquids 867

*<sup>S</sup>*(*λ*)=(<sup>1</sup> <sup>−</sup> *<sup>η</sup>*)2*λ*<sup>3</sup> <sup>+</sup> <sup>6</sup>*η*(<sup>1</sup> <sup>−</sup> *<sup>η</sup>*)*λ*<sup>2</sup> <sup>+</sup> <sup>18</sup>*η*2*<sup>λ</sup>* <sup>−</sup> <sup>12</sup>*η*(<sup>1</sup> <sup>+</sup> <sup>2</sup>*η*).

By expressing Γ as a function of *β*, Henderson *et al*. (3) managed to write the free energy *F* of

∞ ∑ *n*=1

*<sup>ε</sup>* versus *<sup>ρ</sup>*<sup>∗</sup> <sup>=</sup> *ρσ*<sup>3</sup> <sup>=</sup> <sup>6</sup>*<sup>η</sup>*

,

1 + 4*λψ* + 6*λ*2*ψ*<sup>2</sup>

1 + 5*λψ* + 11*λ*2*ψ*<sup>2</sup> + 11*λ*3*ψ*<sup>3</sup>

,

and (1 − *α*<sup>1</sup> + *α*0*ψ*)*φ*<sup>0</sup> = 1.

,

of state, for values of *λ* = 1.8, 3 and 4, are displayed in figure (8b), together with simulation

perturbation theory is noticeably greater than that predicted by the simulation. The reason is that the HCY equation of state involves a truncated series. But, also, we can speculate that the

,

*<sup>λ</sup>*5*φ*<sup>0</sup>

*<sup>λ</sup>*2(<sup>1</sup> <sup>−</sup> *<sup>η</sup>*)

*<sup>λ</sup>*7*φ*<sup>0</sup>

data. A rapid glance at figure (8b) indicates that the critical temperature *T*∗

*vn <sup>n</sup>* (*βε*)

*<sup>n</sup>* , (57)

*<sup>π</sup>* ) predicted by this equation

*<sup>C</sup>* predicted by the

2

where *FHS* is the free energy of the hard-sphere system calculated, for example, with equation

At this stage, we see that the free energy of the HCY potential can be approximated for all values of the triplet (*η*,*ε*, *λ*), using the analytical expressions of the five coefficients *vn* given in footnote (11). From equation (57), we can also deduce the internal energy (Eq. 9), excess pressure (Eq. 12), chemical potential (Eq. 19) and entropy (Eq. 20) of the HCY system. The HCY equation of state, which is based on the perturbation theory and expressed in terms of the relevant features of the potential, is a very handy tool for investigating the thermodynamics of systems governed by an effective hard-sphere interaction plus an attractive tail. Alternatively, the HCY system could be used vicariously to approximate any

*<sup>φ</sup>*<sup>0</sup> <sup>=</sup> *<sup>L</sup>*(*λ*) exp(−*λ*) + *<sup>S</sup>*(*λ*)

*L*(*λ*) = 12*η*[1 + 2*η* + (1 + *η*/2)*λ*],

the system according to the series expansion in powers of *β*:

available interatomic potential for real fluids. The reduced phase diagrams (*T*<sup>∗</sup> = *kBT*

<sup>11</sup> The first terms of the series expansion are:

*<sup>α</sup>*<sup>0</sup> <sup>=</sup> *<sup>L</sup>*(*λ*)

with:

*v*<sup>0</sup> = 0, *<sup>v</sup>*<sup>1</sup> <sup>=</sup> <sup>2</sup>*α*<sup>0</sup> *φ*0 ,

*<sup>v</sup>*<sup>2</sup> <sup>=</sup> <sup>2</sup>*w*(<sup>1</sup> <sup>−</sup> *<sup>α</sup>*<sup>1</sup> <sup>+</sup> *<sup>α</sup>*0*ψ*) *λφ*<sup>0</sup>

*<sup>v</sup>*<sup>4</sup> <sup>=</sup> <sup>4</sup>*w*3(<sup>1</sup> <sup>−</sup> *<sup>α</sup>*<sup>1</sup> <sup>+</sup> *<sup>α</sup>*0*ψ*)

*<sup>v</sup>*<sup>5</sup> <sup>=</sup> <sup>10</sup>*w*4(<sup>1</sup> <sup>−</sup> *<sup>α</sup>*<sup>1</sup> <sup>+</sup> *<sup>α</sup>*0*ψ*)

*<sup>λ</sup>*2(<sup>1</sup> <sup>−</sup> *<sup>η</sup>*)<sup>2</sup> ; *<sup>α</sup>*<sup>1</sup> <sup>=</sup> <sup>12</sup>*η*(<sup>1</sup> <sup>+</sup> *<sup>λ</sup>*/2)

*<sup>v</sup>*<sup>3</sup> <sup>=</sup> <sup>2</sup>*w*2(<sup>1</sup> <sup>−</sup> *<sup>α</sup>*<sup>1</sup> <sup>+</sup> *<sup>α</sup>*0*ψ*)(<sup>1</sup> <sup>+</sup> <sup>3</sup>*λψ*) *<sup>λ</sup>*3*φ*<sup>0</sup>

*<sup>λ</sup>*3(<sup>1</sup> <sup>−</sup> *<sup>η</sup>*)<sup>2</sup> ,

*<sup>F</sup>* <sup>=</sup> *FHS* <sup>−</sup> *NkBT*

(25), and where the first five terms *vn* have expressions derived by the authors11.

with:

#### **5.6 Application to the Yukawa attractive potential**

As an application of the perturbation theory, consider the hard-core attractive Yukawa potential (HCY). This potential consists of the hard-sphere potential *uHS*(*d*) and the attractive perturbation *u*1(*r*). Its traditional analytical form is:

$$\mu(r) = \begin{cases} \infty & r < d \\ -\left(\varepsilon d / r\right) \exp\left[-\lambda \left(r/d - 1\right)\right] r \ge d \end{cases} \tag{55}$$

where *d* is the diameter of hard spheres, *ε* the potential value at *r* = *d* and *λ* a parameter that measures the rate of exponential decay. The shape of this potential is portrayed in figure (8a)

Fig. 8. (a) Representation of the HCY potential for three values of *λ* (1.8, 3.0, 4.0). (b) Phase diagram of the HCY potential for the same three values of *λ* (1.8, 3.0, 4.0).

for three different values of *λ*. For guidance, note that a value of *λ* ∼ 1.8 predicts the thermodynamic properties of liquid rare gases quite so well as the Lennard-Jones potential. By contrast, a value of *λ* ∼ 8 enables us to obtain the thermodynamic properties of colloidal suspensions or globular proteins. This wide possible range of *λ* values explains why the HCY potential has been used in many applications and has been the subject of numerous theoretical studies after its analytical solution was obtained by Waisman (2). Without going into details of the resolution of the problem that involves the microscopic theory of liquids, indicate that the solution reduces to determining a fundamental parameter Γ as the root of the quartic equation (23):

$$
\Gamma(1+\lambda\Gamma)(1+\psi\Gamma)^2+\beta\epsilon w=0,\tag{56}
$$

where *ε* and *λ* are the two parameters of the HCY potential, and *w* and *ψ* two additional parameters explicitly depending on *λ* and *η* by the following relations:

$$\begin{aligned} w &= \frac{6\eta}{\phi\_0^2}, \\ \psi &= \lambda^2 (1-\eta)^2 \frac{[1-\exp(-\lambda)]}{L(\lambda)\exp(-\lambda) + S(\lambda)} \\ &- 12\eta (1-\eta) \frac{[1-\lambda/2 - (1+\lambda/2)\exp(-\lambda)]}{L(\lambda)\exp(-\lambda) + S(\lambda)} \end{aligned}$$

with:

28 Thermodynamics book 1

As an application of the perturbation theory, consider the hard-core attractive Yukawa potential (HCY). This potential consists of the hard-sphere potential *uHS*(*d*) and the attractive

−(*εd*/*r*) exp[−*λ*(*r*/*d* − 1)]

where *d* is the diameter of hard spheres, *ε* the potential value at *r* = *d* and *λ* a parameter that measures the rate of exponential decay. The shape of this potential is portrayed in figure (8a)

Fig. 8. (a) Representation of the HCY potential for three values of *λ* (1.8, 3.0, 4.0). (b) Phase

for three different values of *λ*. For guidance, note that a value of *λ* ∼ 1.8 predicts the thermodynamic properties of liquid rare gases quite so well as the Lennard-Jones potential. By contrast, a value of *λ* ∼ 8 enables us to obtain the thermodynamic properties of colloidal suspensions or globular proteins. This wide possible range of *λ* values explains why the HCY potential has been used in many applications and has been the subject of numerous theoretical studies after its analytical solution was obtained by Waisman (2). Without going into details of the resolution of the problem that involves the microscopic theory of liquids, indicate that the solution reduces to determining a fundamental parameter Γ as the root of the quartic equation

where *ε* and *λ* are the two parameters of the HCY potential, and *w* and *ψ* two additional

*L*(*λ*) exp(−*λ*) + *S*(*λ*)

[1 − *λ*/2 − (1 + *λ*/2) exp(−*λ*)] *<sup>L</sup>*(*λ*) exp(−*λ*) + *<sup>S</sup>*(*λ*) ,

diagram of the HCY potential for the same three values of *λ* (1.8, 3.0, 4.0).

parameters explicitly depending on *λ* and *η* by the following relations:

−12*η*(1 − *η*)

*<sup>ψ</sup>* <sup>=</sup> *<sup>λ</sup>*2(<sup>1</sup> <sup>−</sup> *<sup>η</sup>*)<sup>2</sup> [<sup>1</sup> <sup>−</sup> exp(−*λ*)]

*<sup>w</sup>* <sup>=</sup> <sup>6</sup>*<sup>η</sup> φ*2 0 ,

(23):

*r* < *d*

Γ(1 + *λ*Γ)(1 + *ψ*Γ)<sup>2</sup> + *βεw* = 0, (56)

*<sup>r</sup>* <sup>≥</sup> *<sup>d</sup>*, (55)

∞

**5.6 Application to the Yukawa attractive potential**

perturbation *u*1(*r*). Its traditional analytical form is:

*u*(*r*) =

$$\begin{aligned} \phi\_0 &= \frac{L(\lambda) \exp(-\lambda) + S(\lambda)}{\lambda^3 (1 - \eta)^2}, \\ L(\lambda) &= 12\eta [1 + 2\eta + (1 + \eta/2)\lambda], \\ S(\lambda) &= (1 - \eta)^2 \lambda^3 + 6\eta (1 - \eta) \lambda^2 + 18\eta^2 \lambda - 12\eta (1 + 2\eta). \end{aligned}$$

By expressing Γ as a function of *β*, Henderson *et al*. (3) managed to write the free energy *F* of the system according to the series expansion in powers of *β*:

$$F = F\_{HS} - \frac{Nk\_B T}{2} \sum\_{n=1}^{\infty} \frac{v\_n}{n} \left(\beta \varepsilon\right)^n,\tag{57}$$

where *FHS* is the free energy of the hard-sphere system calculated, for example, with equation (25), and where the first five terms *vn* have expressions derived by the authors11.

At this stage, we see that the free energy of the HCY potential can be approximated for all values of the triplet (*η*,*ε*, *λ*), using the analytical expressions of the five coefficients *vn* given in footnote (11). From equation (57), we can also deduce the internal energy (Eq. 9), excess pressure (Eq. 12), chemical potential (Eq. 19) and entropy (Eq. 20) of the HCY system.

The HCY equation of state, which is based on the perturbation theory and expressed in terms of the relevant features of the potential, is a very handy tool for investigating the thermodynamics of systems governed by an effective hard-sphere interaction plus an attractive tail. Alternatively, the HCY system could be used vicariously to approximate any available interatomic potential for real fluids.

The reduced phase diagrams (*T*<sup>∗</sup> = *kBT <sup>ε</sup>* versus *<sup>ρ</sup>*<sup>∗</sup> <sup>=</sup> *ρσ*<sup>3</sup> <sup>=</sup> <sup>6</sup>*<sup>η</sup> <sup>π</sup>* ) predicted by this equation of state, for values of *λ* = 1.8, 3 and 4, are displayed in figure (8b), together with simulation data. A rapid glance at figure (8b) indicates that the critical temperature *T*∗ *<sup>C</sup>* predicted by the perturbation theory is noticeably greater than that predicted by the simulation. The reason is that the HCY equation of state involves a truncated series. But, also, we can speculate that the

<sup>11</sup> The first terms of the series expansion are:

$$\begin{aligned} v\_0 &= 0, \\ v\_1 &= \frac{2\alpha\_0}{\Phi\_0}, \\ v\_2 &= \frac{2w(1 - \mathfrak{a}\_1 + \mathfrak{a}\_0\mathfrak{\psi})}{\lambda\Phi\_0}, \\ v\_3 &= \frac{2w^2(1 - \mathfrak{a}\_1 + \mathfrak{a}\_0\mathfrak{\psi})\left(1 + 3\lambda\mathfrak{\psi}\right)}{\lambda^3\Phi\_0}, \\ v\_4 &= \frac{4w^3(1 - \mathfrak{a}\_1 + \mathfrak{a}\_0\mathfrak{\psi})\left(1 + 4\lambda\mathfrak{\psi} + 6\lambda^2\mathfrak{\psi}^2\right)}{\lambda^5\Phi\_0}, \\ v\_5 &= \frac{10w^4(1 - \mathfrak{a}\_1 + \mathfrak{a}\_0\mathfrak{\psi})\left(1 + 5\lambda\mathfrak{\psi} + 11\lambda^2\mathfrak{\psi}^2 + 11\lambda^3\mathfrak{\psi}^3\right)}{\lambda^7\Phi\_0}, \end{aligned}$$

with:

$$\mathfrak{a}\_0 = \frac{L(\lambda)}{\lambda^2 (1 - \eta)^2} \quad ; \quad \mathfrak{a}\_1 = \frac{12\eta (1 + \lambda/2)}{\lambda^2 (1 - \eta)} \quad \text{and} \quad (1 - \mathfrak{a}\_1 + \mathfrak{a}\_0 \mathfrak{y}) \mathfrak{d}\_0 = 1.$$

of Simple Liquids 31

Thermodynamic Perturbation Theory of Simple Liquids 869

expressions for the hard-sphere free energy *FHS*(*d*1, *d*2) and subsequent thermodynamic quantities are readily generalized for a mixture of hard spheres of different diameters *d*<sup>1</sup> and *d*2. Therefore, if one estimates that a real binary mixture behaves like binary hard-sphere fluid, the variational method consists of minimizing *FHS*(*d*1, *d*2) versus the two hard-sphere diameters, and by taking the resulting minimum upper bound as an approximation to the free energy. The variational method is not limited neither to hard-sphere reference system nor to a specific fluid. It should be stressed that systems like pure liquid metals composed of ions embedded in an electron gas can not be treated as a binary mixture, in the strict sense of the word, but must be suitably reduced to a one-component system of pseudoions. In contrast, a liquid metal made up of two species of pseudoions forms a binary mixture, for which the approach originally developed for simple liquids can be applied. What is it that determines whether or not two metals will mix to form an alloy is a crucial issue to be answered by

I would like to express my acknowledgement to Jean-Marc Bomont for its stimulating

[3] D. Henderson, G. Stell and E. Waisman, J. Chem. Phys. 62, 4247 (1975). D. Henderson, L.

[4] See, for instance, the book of J. P. Hansen and I. R. McDonald, *Theory of Simple Liquids*,

[6] E. Chacon, M. Reinaldo-Falagan, E. Velasco and P. Tarazona, Phys. Rev. Lett. 87, 166101 (2001). D. Li and S. A. Rice, J. Phys. Chem. B. 108, 19640 (2004). J. M. Bomont and J. L.

[8] J. G. Kirkwood, J. Chem. Phys. 3, 300 (1935). J. G. Kirkwood and E. Monroe, J. Chem. Phys. 9, 514 (1941). See also D. A. McQuarrie, *Statistical Mechanics*, Harper and Row,

[9] R. E. Nettleton and M. S. Green, J. Chem. Phys. 29, 1365 (1958). H. J. Raveché, J. Chem.

[10] B. J. Alder and T. E. Wainwright, Phys. Rev. 127, 359 (1962); F. H. Ree and W. G. Hoover,

[20] G. A. Mansoori and F. B. Canfield, J. Chem. Phys. 51, 4967 (1969); 51 5295 (1969); 53, 1618

[21] J. D. Weeks, D. Chandler and H. C. Andersen, J. Chem. Phys. 54, 5237 (1971); 55, 5422

[5] W. G. Hoover and F. H. Ree, J. Chem. Phys. 47, 4873 (1967); 49, 3609 (1968).

[7] N. Jakse and J. L. Bretonnet, J. Phys.: Condens. Matter 15, S3455 (2003).

[11] N. F. Carnahan and K. E. Starling, J. Chem. Phys. 51, 635 (1969).

[17] J. A. Barker and D. Henderson, J. Chem. Phys. 47, 2856 (1967). [18] G. A. Mansoori and F. B. Canfield, J. Chem. Phys. 51, 4958 (1969). [19] J. A. Barker and D. Henderson, Ann. Rev. Phys. Chem. 23, 439 (1972).

[14] N. W. Ashcroft and J. Lekner, Phys. Rev. 156, 83 (1966).

[12] D. A. McQuarrie, *Statistical Mechanics*, Harper and Row, New York (1976). [13] A. J. Greenfield, J. Wellendorf and N. Wiser, Phys. Rev. A 4, 1607 (1971).

[1] M. H. J. Hagen and D. Frenkel, J. Chem. Phys. 101, 4093 (1994).

Blum and J. P. Noworyta, J. Chem. Phys. 102, 4973 (1995).

Bretonnet, J. Chem. Phys. 124, 054504 (2006).

thermodynamic perturbation theory (27).

[2] E. Waisman, Mol. Phys. 25, 45 (1973).

Acad. Press., London (2006).

New-York (1976), p. 263.

J. Chem. Phys. 40, 939 (1964).

[15] L. Verlet, Phys. Rev. 165, 201 (1968). [16] R. Zwanzig, J. Chem. Phys. 22, 1420 (1954).

Phys. 55, 2242 (1971).

(1970).

(1971).

discussions.

**7. References**

correlation functions, inherent in the calculations of the pressure and the chemical potential, do not represent correctly the growing correlation lengths in approaching the critical point. Looking at the evolution of the binodal lines as a function of the rate of decay *λ* of the Yukawa potential, we remark that higher the critical temperature is, lower *λ* is. At the same time, the domain below the binodal line shrinks. In contrast, when *λ* increases (i. e., the attraction range of the potential becomes shorter), the critical temperature decreases, and the liquid and gaseous phases become indistinguishable. For the hard-sphere potential, as a border case (*λ* → ∞), there is no longer gas-liquid phase transition. On the other hand, one can see that the phase diagrams obtained with the perturbation theory agree with simulation data more favorably for the vapor branch than for the liquid branch. This is not surprising in the extent that the perturbation theory works better for low densities. Lastly, it should be mentioned that the structure and thermodynamic properties of the HCY potential have been extensively studied for the two last decades, as much by computer simulations and integral equation theory as by means of perturbation theory. Note that many other studies of the HCY potential with these various methods are available in literature (24).
