Thermodynamic Perturbation Theory of Simple Liquids

This chapter is an introduction to the thermodynamics of systems, based on the correlation function formalism, which has been established to determine the thermodynamic properties of simple liquids. The article begins with a preamble describing few general aspects of the liquid state, among others the connection between the phase diagram and the pair potential u(r), on one hand, and between the structure and the pair correlation function g(r), on the other hand. The pair correlation function is of major importance in the theory of liquids at equilibrium, because it is required for performing the calculation of the thermodynamic properties of systems modeled by a given pair potential. Then, the article is devoted to the expressions useful for calculating the thermodynamic properties of liquids, in relation with the most relevant features of the potential, and provides a presentation of the perturbation theory developed in the four last decades. The thermodynamic perturbation theory is founded on a judicious separation of the pair potential into two parts. Specifically, one of the greatest successes of the microscopic theory has been the recognition of the quite distinct roles played by the repulsive and attractive parts of the pair potential in predicting many properties of liquids. Much attention has been paid to the hard-sphere potential, which has proved very efficient as natural reference system because it describes fairly well the local order in liquids. As an example, the Yukawa attractive potential is also mentioned.


Introduction
This chapter is an introduction to the thermodynamics of systems, based on the correlation function formalism, which has been established to determine the thermodynamic properties of simple liquids. The article begins with a preamble describing few general aspects of the liquid state, among others the connection between the phase diagram and the pair potential u(r), on one hand, and between the structure and the pair correlation function g(r),o nt h e other hand. The pair correlation function is of major importance in the theory of liquids at equilibrium, because it is required for performing the calculation of the thermodynamic properties of systems modeled by a given pair potential. Then, the article is devoted to the expressions useful for calculating the thermodynamic properties of liquids, in relation with the most relevant features of the potential, and provides a presentation of the perturbation theory developed in the four last decades. The thermodynamic perturbation theory is founded on a judicious separation of the pair potential into two parts. Specifically, one of the greatest successes of the microscopic theory has been the recognition of the quite distinct roles played by the repulsive and attractive parts of the pair potential in predicting many properties of liquids. Much attention has been paid to the hard-sphere potential, which has proved very efficient as natural reference system because it describes fairly well the local order in liquids. As an example, the Yukawa attractive potential is also mentioned.

An elementary survey 2.1 The liquid state
The ability of the liquids to form a free surface differs from that of the gases, which occupy the entire volume available and have diffusion coefficients (∼ 0, 5 cm 2 s −1 ) of several orders of magnitude higher than those of liquids (∼ 10 −5 cm 2 s −1 ) or solids (∼ 10 −9 cm 2 s −1 ). Moreover, if the dynamic viscosity of liquids (between 10 −5 Pa.s and 1 Pa.s) is so lower compared to that of solids, it is explained in terms of competition between configurational and kinetic processes. Indeed, in a solid, the displacements of atoms occur only after the breaking of the bonds that keep them in a stable configuration. At the opposite, in a gas, molecular transport is a purely kinetic process perfectly described in terms of exchanges of energy and momentum. In a liquid, the continuous rearrangement of particles and the molecular transport combine together in appropriate proportion, meaning that the liquid is an intermediate state between the gaseous and solid states. The characterization of the three states of matter can be done in an advantageous manner by comparing the kinetic energy and potential energy as it is done in figure (1). The nature and intensity of forces acting between particles are such that the particles tend to attract each other at great distances, while they repel at the short distances. The particles are in equilibrium when the attraction and repulsion forces balance each other. In gases, the kinetic energy of particles, whose the distribution is given by the Maxwell velocity distribution, is located in the region of unbound states. The particles move freely on trajectories suddenly modified by binary collisions; thus the movement of particles in the gases is essentially an individual movement. In solids, the energy distribution is confined within the potential well. It follows that the particles are in tight bound states and describe harmonic motions around their equilibrium positions; therefore the movement of particles in the solids is essentially a collective movement. When the temperature increases, the energy distribution moves towards high energies and the particles are subjected to anharmonic movements that intensify progressively. In liquids, the energy distribution is almost entirely located in the region of bound states, and the movements of the particles are strongly anharmonic. On approaching the critical point, the energy distribution shifts towards the region of unbound states. This results in important fluctuations in concentrations, accompanied by the destruction and formation of aggregates of particles. Therefore, the movement of particles in liquids is thus the result of a combination of individual and collective movements. When a crystalline solid melts, the long-range order of the crystal is destroyed, but a residual local order persits on distances greater than several molecular diameters. This local order into liquid state is described in terms of the pair correlation function, g(r)= ρ(r) ρ ∞ , which is defined as the ratio of the mean molecular density ρ(r), at a distance r from an arbitrary molecule, to the bulk density ρ ∞ .I fg(r) is equal to unity everywhere, the fluid is completely disordered, like in diluted gases. The deviation of g(r) from unity is a measure of the local order in the arrangement of near-neighbors. The representative curve of g(r) for a liquid is formed of maxima and minima rapidly damped around unity, where the first maximum corresponds 4 Thermodynamics book 1 of statistical thermodynamics providing a link between the microscopic description of liquids and classical thermodynamic functions. Then, it will be given an account of the thermodynamic perturbation theory with the analytical expressions required for calculating the thermodynamic properties. Finally, the HCY model, which is founded on the perturbation theory, will be presented in greater detail for investigating the thermodynamics of liquids. Thus, a review of the thermodynamic perturbation theory will be set up, with a special effort towards the pedagogical aspect. We hope that this paper will help readers to develop their inductive and synthetic capacities, and to enhance their scientific ability in the field of thermodynamic of liquids. It goes without saying that the intention of the present paper is just to initiate the readers to that matter, which is developed in many standard textbooks (4).

Phase stability limits versus pair potential
One success of the numerical simulation was to establish a relationship between the shape of the pair potential and the phase stability limits, thus clarifying the circumstances of the liquid-solid and liquid-vapor phase transitions. It has been shown, in particular, that the hard-sphere (HS) potential is able to correctly describe the atomic structure of liquids and predict the liquid-solid phase transition (5). By contrast, the HS potential is unable to describe the liquid-vapor phase transition, which is essentially due to the presence of attractive forces of dispersion. More specifically, the simulation results have shown that the liquid-solid phase transition depends on the steric hindrance of the atoms and that the coexistence curve of liquid-solid phases is governed by the details of the repulsive part of potential. In fact, this was already contained in the phenomenological theories of melting, like the Lindemann theory that predicts the melting of a solid when the mean displacement of atoms from their equilibrium positions on the network exceeds the atomic diameter of 10%. In other words, a substance melts when its volume exceeds the volume at0Kof30%. In restricting the discussion to simple centrosymmetric interactions from the outset, it is necessary to consider a realistic pair potential adequate for testing the phase stability limits. The most natural prototype potential is the Lennard-Jones (LJ) potential given by where the parameters m and n are usually taken to be equal to 12 and 6, respectively. Such a functional form gives a reasonable representation of the interactions operating in real fluids, where the well depth ε LJ and the collision diameter σ LJ are independent of density and temperature. Figure (2a) displays the general shape of the Lennard-Jones potential (m − n) corresponding to equation (1). Each substance has its own values of ε LJ and σ LJ so that, in reduced form, the LJ potentials have not only the same shape for all simple fluids, but superimpose each other rigorously. This is the condition for substances to conform to the law of corresponding states. Figure (2b) represents the diagram p(T) of a pure substance. We can see how the slope of the coexistence curve of solid-liquid phases varies with the repulsive part of potential: the higher the value of m, the steeper the repulsive part of the potential (Fig. 2a) and, consequently, the more the coexistence curve of solid-liquid phases is tilted (Fig. 2b).
We can also remark that the LJ potential predicts the liquid-vapor coexistence curve, which begins at the triple point T and ends at the critical point C. A detailed analysis shows that the length of the branch TC is proportional to the depth ε of the potential well. As an example, for rare gases, it is verified that (T C − T T )k B ≃ 0, 55 ε. It follows immediately from this condition that the liquid-vapor coexistence curve disappears when the potential well is absent (ε = 0).  The value of the slope of the branch TC also depends on the attractive part of the potential as shown by the Clausius-Clapeyron equation: where L vap is the latent heat of vaporization at the corresponding temperature T vap and (V vap − V liq ) is the difference of specific volumes between vapor and liquid. To evaluate the slope dp dT of the branch TC at ambient pressure, we can estimate the ratio L vap T vap with Trouton's rule ( L vap T vap ≃ 85 J.K −1 .mol −1 ), and the difference in volume (V vap − V liq ) in terms of width of the potential well. Indeed, in noting that the quantity (V vap − V liq ) is an increasing function of the width of potential well, which itself increases when n decreases, we see that, for a given well depth ε, the slope of the liquid-vapor coexistence curve decreases as n decreases. For liquid metals, it should be mentioned that the repulsive part of the potential is softer than for liquid rare gases. Moreover, even if ε is slightly lower for metals than for rare gases, the is much higher (between 2 and 4), which explains the elongation of the TC curve compared to that of rare gases. It is worth also to indicate that some flat-bottomed potentials (6) are likely to give a good description of the physical properties of substances that have a low value of the ratio T T T C . Such a potential is obviously not suitable for liquid rare gases, whose ratio T T T C ≃ 0, 56, or for organic and inorganic liquids, for which 0, 25 < T T T C < 0, 45. In return, it might be useful as empirical potential for metals with low melting point such as mercury, gallium, indium, tin, etc., the ratio of which being T T T C < 0, 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 (λ ∼ 0, 1 nm). At this radiation corresponds a photon energy (hν = hc λ ∼ 10 4 eV), much larger than the mean energy of atoms that is of the order of few k B T, namely about 10 −1 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 σ = 4πb 2 , where b is a parameter approximately equal to the radius of the core (∼ 10 −14 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.

Structure factor and pair correlation function
When a photon of wave vector k = 2π λ u interacts with an atom, it is deflected by an angle θ and the wave vector of the scattered photon is k ′ = 2π λ u ′ , where u and u ′ are unit vectors. If the scattering is elastic it results that |k ′ | = |k| , because E ∝ k 2 = cte, and that the scattering vector (or transfer vector) q is defined by the Bragg law: 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: 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: 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 P N (r N , p N ) of the canonical ensemble 1 . Therefore, the statistical average can be written by using the distribution 1 It seems useful to remember that the probability density function in the canonical ensemble is: where H N (r N , p N )=∑ p 2 2m + U(r N ) is the Hamiltonian of the system, β = 1 k B T and Q N (V, T) the partition function defined as: , 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 Z N (V, T), which is expressed in terms of the total potential energy U(r N ). They read: Besides, the partition function Q N (V, T) allows us to determine the free energy F according to the relation: The reader is advised to consult statistical-physics textbooks for further details.
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 in the statistical average as 3 : 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): 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 function 4 , 2 It should be stressed that the distribution function ρ (2) N r 2 is expressed as: 3 To evaluate an integral of the form: one must use the spherical coordinates by placing the vector q along the z axis, where θ =(q, r). Thus, the integral reads: with μ = cos θ and dμ = − sin θdθ. It follows that: qr g(r)r 2 dr. 4 The generalization of the Fourier transform of the Dirac delta function to three dimensions is: and the inverse transform is: 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: 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: 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: 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)=ρk B Tχ T .

847
Thermodynamic Perturbation Theory of Simple Liquids www.intechopen.com

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 : where the partition function Q N (V, T) depends on the configuration integral Z N (V, T) and on the thermal wavelength Λ, in accordance with the equations given in footnote (1). The derivative of ln Q N (V, T) with respect to T can be written: 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 )=∑ i ∑ j>i u(r ij ). 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, permitting us to write: where the mean value is expressed in terms of the pair correlation function as: For a homogeneous and isotropic fluid, one can perform the change of variables R = r 1 and r = r 1 − r 2 , where R and r describe the system volume, and write the expression of internal energy in the integral form: 848 Thermodynamics -Interaction Studies -Solids, Liquids and Gases www.intechopen.com Thermodynamic Perturbation Theory of Simple Liquids 11 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.

Pressure
The expression of the pressure is obtained in the same way that the internal energy, in considering the equation: The derivation of the configuration integral with respect to volume requires using the reduced variable X = r V that allows us to find the dependence of the potential energy U(r N ) versus volume. Indeed, if the volume element is dr = VdX, the scalar variable dr = V 1/3 dX leads to the derivative: In view of this, the configuration integral and its derivative with respect to V are written in the following forms with reduced variables: 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: so that the expression of the pressure becomes: 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:

849
Thermodynamic Perturbation Theory of Simple Liquids where the mean value is expressed with the pair correlation function by: For a homogeneous and isotropic fluid, one can perform the change of variables R = r 1 and r = r 1 − r 2 , and simplify the expression of pressure as: 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 ).

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: From footnote (1), the free energy F is written: so that the chemical potential can be simplified as: 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: 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 850 Thermodynamics -Interaction Studies -Solids, Liquids and Gases www.intechopen.com Thermodynamic Perturbation Theory of Simple Liquids 13 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: Under these conditions, the configuration integrals for a total coupling (λ = 1) and a zero coupling (λ = 0) are respectively: These expressions are then used to calculate the logarithm of the ratio of configuration integrals in equation (13): But with the configuration integral Z N (V, T, λ), in which potential energy is given by equation (14), we can easily evaluate the partial derivatives ∂Z N ∂λ and ∂ ln Z N ∂λ . In particular, with the result of the footnote (1), we can write ∂ ln Z N ∂λ as a function of the pair correlation function as: In addition, if the fluid is homogeneous and isotropic, the above relation simplifies under the following form: 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: 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): 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).

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: 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 is a step function (Fig. 4) whose derivative with respect to r is the Dirac delta function, i. e.: In substituting ∂u ∂r , taken from the previous relation, in equation (12) we find the expression of the pressure: It is important to recall that, for moderately dense gases, the pressure is usually expressed under the form of the virial expansion The first term of the last equality represents the contribution of the ideal gas, and the excess pressure p ex comes from the interactions between particles. They are written: 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: Note that the first 6 coefficients of the excess pressure p ex 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: In combining the first and second derivatives of the geometric series ∑ ∞ k=1 η k , it is found that equation (23) can be transformed into a rational fraction 5 enabling the deduction of the excess pressure in the form: Consequently, the equation of state of the hard-sphere fluid is written with excellent precision as: 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 = 3 2 Nk B T. As for the free energy F, it is determined by integrating the pressure over volume with the equation: where F GP is the free energy of ideal gas (cf. footnote 1, with Z N (V, T)=V N ): and F ex the excess free energy, calculated by integrating equation (23) as follows: 5 To obtain the rational fraction, one must decompose the sum as: and combine the geometric series ∑ ∞ k=1 η k with its first and second derivatives: to see appear the relation: But, with equation (22) that gives dV dη = − V η , F ex is then reduced to the series expansion: Like the pressure, this expansion is written as a rational function by combining the geometric series ∑ ∞ k=1 η k with its first derivative 6 . The expression of the excess free energy is: 2 , and the free energy of the hard-sphere fluid reduces to the following form: Now, the entropy is obtained using the same method of calculation, by deriving the free energy with respect to temperature: where S GP is the entropy of the ideal gas given by the Sackur-Tetrode equation: and where the excess entropy S ex arises from the relation: hence the expression of the entropy of the hard-sphere fluid: Finally, combining equations (25) and (24), with the help of equation (20), one reaches the chemical potential of the hard-sphere fluid that reads: 6 Indeed, the identity: is yet written:

855
Thermodynamic Perturbation Theory of Simple Liquids

www.intechopen.com
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, g HS (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.

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.  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 S HS (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 methods 7 . 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 g HS (r) of the hard-sphere system and not that of the real system.

Zwanzig method
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 ),i st h e 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: 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: By replacing the potential energy U(r N ) in the expression of the configuration integral (cf. footnote 1), one gets: The first integral represents the configuration integral Z (0) N (V, T) of the reference system, and the remaining term can be regarded as the average value of the quantity exp −βλU 1 (r N ) ,so that the previous relation can be put under the general form: where ... 0 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 7 The interested reader will find all useful adjuncts in the books either by J. P. Hansen 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: Since the perturbation U 1 (r N ) is small, exp (−βλU 1 ) can be expanded in series, so that the statistical average exp −βλU 1 (r N ) 0 , calculated on the reference system, is expressed as: 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 −βλU 1 (r N ) 0 in Taylor series. After simplifications, equation (31) reduces to:

Now if we set:
we find that: ln exp −βλU 1 (r N ) The contribution of the perturbation (Eq. 31) to the free energy is then written in the compact form: and the expression of the free energy F of the real system is found by substituting equation (36) into equation (30), as follows: 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.

858
Thermodynamics -Interaction Studies -Solids, Liquids and Gases www.intechopen.com Thermodynamic Perturbation Theory of Simple Liquids 21

Van der Waals equation
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 1 is to admit that the total potential energy may be decomposed into a sum of pair potentials in the form: Therefore, the free energy of the perturbation to zero order in β is given by equation (33), that is to say: 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 , the expression of c 1 is simplified as: The integral in between the braces is then expressed as a function of the pair correlation function (cf. footnote 2), and c 1 reduces to: 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 particle 8 is b = 2 3 πσ 3 and the unoccupied volume is (V − Nb). With these simplifications in mind, the configuration integral and free energy of the reference system are respectively (cf. footnote 1): 8 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 4 3 πσ 3 must be assigned to each particle (Fig.  6b).  As for the coefficient c 1 (Eq. 38), it is simplified as: Therefore, the expression of free energy (Eq. 37) corresponding to the model of van der Waals is: and the van der Waals equation of state reduces to: 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 g HS (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.

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 ) , ..., (r i , r i+1 ),..., and assimilate the perturbating elemental potential in each interval by a constant. Assuming that the perturbating potential in the interval (r i , r i+1 ) is u 1 (r i ) and that the number of atoms subjected to this potential is N i , the total perturbation can be written as the sum of elemental potentials: before substituting it in the configuration integral. The advantage of this method is to calculate the coefficients c n 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 (r i ) is constant in the interval which it belongs, so that we can write U 1 2 0 = ∑ i ∑ j N i 0 N j 0 u 1 (r i )u 1 (r j ).I n view of this, the coefficient c 2 defined by equation (34) is: With the local compressibility approximation (LC), where ρ and g 0 depend on p, the expression of c 2 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 3 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: 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, 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: By comparing equations (37) and (42), we see that the first derivative ∂F ∂λ coincides with c 1 and the second derivative 1 2 ∂ 2 F ∂λ 2 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 ∂F ∂α , whose the explicit form given without proof, reads: Since the Barker and Henderson prescription is based on the proposal to cancel the term ∂F ∂α , 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 of F in α 2 and αλ are negligible compared to the term in λ 2 .

862
Thermodynamics -Interaction Studies -Solids, Liquids and Gases www.intechopen.com Thermodynamic Perturbation Theory of Simple Liquids 25 Therefore, using equation (38) to evaluate c 1 and equation (40) to evaluate c 2 , the expression of the free energy F of the real system (Eq. 42) is: where the first term on the RHS represents the free energy of the hard-sphere system (Eq. 25), and the partial derivative Recall that the pair correlation function of hard-sphere system, g HS (η; r), must be only determined numerically. It depends on the density ρ and diameter d via the packing fraction η(= π 6 ρd 3 ). Since g HS (η; r)=0 when r < d, either 0 or d can be used as lower limit of integration in equation (44).

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 1 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 0 as: After replacing U 1 in equation (30), we obtain: However, considering the series expansion of an exponential, the above relation is transformed into the so-called Gibbs-Bogoliubov inequality 9 : 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 η = π 6 ρd 3 ) is considered as the best estimate of the free energy of the real system. Its expression is: 9 The Mac-Laurin series of the exponential naturally leads to the inequality: For n = 0, the last term of equation (45) behaves as: since the mean value of the deviation, ΔU 1 0 , is zero. where F 0 = F HS is given by equation (25) and U 1 0 = c 1 by equation (38). Since g HS (η; r)=0 when r < d and u 1 (r)=u(r) − u HS (r)=u(r) when r ≥ d, equation (47) can be written as: 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).

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 − ∂u ∂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: 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  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, F HS , 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: where y HS (η; r) is the cavity function of great importance in the microscopic theory of liquids 10 . It is defined by means of the pair correlation function g HS (η; r) and the Boltzmann factor exp [−βu HS (d)] of the hard-sphere potential as: 10 The cavity function y HS (η; r) does not exist in analytical form. It must be calculated at all reduced densities (ρd 3 ) by solving an integro-differential equation. 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 y HS (η; r) for all potentials, and to express the pair correlation functions of each potential as a function of y HS (η; 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: The suggestion of WCA for choosing the hard-sphere diameter d is to cancel ∂F ∂α , which amounts to solving the nonlinear equation given by equation (50), i. e.: 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 g HS (η; r) and g 0 (r) drawn, respectively, from equations (51) and (52). It reads: According to equation (7), the previous expression is equivalent to the equality S 0 (0)= S HS (0), meaning the equality between the isothermal compressibility of the repulsive potential u 0 (r) and that of the hard-sphere potential u HS (d).
Ultimately, the expression of free energy (Eq. 28) of the real system is the sum of F HS , calculated with the value of d issued from equation (53), and the term ∂F ∂λ = c 1 , calculated with equation (38) in using the pair correlation function g 0 (r) ≃ exp [−βu 0 (r)] y HS (η; 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 y HS (η; 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 α 2 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 F HS , c 1 and c 2 are calculated with the same reference system, whereas the WCA treatment is not coherent because it uses the hard-sphere system to calculate F HS 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 2 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).

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 u HS (d) and the attractive perturbation u 1 (r). Its traditional analytical form is: 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) 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): where ε and λ are the two parameters of the HCY potential, and w and ψ two additional parameters explicitly depending on λ and η by the following relations: with: 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 β: where F HS is the free energy of the hard-sphere system calculated, for example, with equation (25), and where the first five terms v n have expressions derived by the authors 11 . 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 v n 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 * = k B T ε versus ρ * = ρσ 3 = 6η π ) 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 * C 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 11 The first terms of the series expansion are: with:

867
Thermodynamic Perturbation Theory of Simple Liquids www.intechopen.com 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).

Concluding remarks
In the foregoing, we started with a brief discussion of the liquid state, compared with the gaseous and solid states. Among other things, we mentioned that one of the most successes of the microscopic theory of liquids has been to emphasize the importance of the pair potential u(r) to describe a wide variety of physical properties of liquids. Then, we learned a rudimentary knowledge of the structure of liquids. The latter is essentially described by the static structure factor S(q), which can be measured directly by elastic scattering of neutrons or X-rays. But S(q), or more precisely the pair correlation function g(r), can also be determined by numerical simulations (known as virtual experiments) and with models of the microscopic theory of liquids once the nature of the pair potential u(r) is known. The comparison between simulation results and those of analytical models essentially allows us to test the models, whereas the comparison of simulation results with experimental results is the ultimate test to judge the efficiency of the pair potential. Subsequently, we described in minute detail the calculations of the thermodynamic properties in terms of the pair potential and pair correlation function. In particular, with the pressure p and the chemical potential μ to apply equilibrium conditions, it has been shown that we are in a position to determine theoretical estimates for the liquid-vapor coexistence curve. Finally, we got down to basics of the thermodynamic perturbation theory, before presenting the liquid-vapor coexistence curve for the HCY fluid. In a 1990 time warp, it became apparent that the usefulness of thermodynamic perturbation theory was on the decline compared to the more powerful simulation methods. This was primarily due to the rapid increase in the power of computers. At the same time, the integral equation theory enjoyed renewed popularity with the extensively employed concept of thermodynamic consistency (25), but this aspect goes beyond the scope of this short review. Incidentally, we recall that the thermodynamic consistency consists to adjusting the isothermal compressibility obtained by two different routes. Nevertheless, the thermodynamic perturbation theory remains undoubtedly the most tractable approach to predict the thermodynamic properties of liquids. Let us just quote few articles containing investigations of the HCY fluid with the integral equation theory (26). As one would expect, the thermodynamic perturbation theory described previously is also important for mixtures. The thermodynamic properties established for studying pure liquids can be applied to binary mixtures, with only few modifications. Specifically, analytic 868 Thermodynamics -Interaction Studies -Solids, Liquids and Gases www.intechopen.com Thermodynamic Perturbation Theory of Simple Liquids 31 expressions for the hard-sphere free energy F HS (d 1 , d 2 ) and subsequent thermodynamic quantities are readily generalized for a mixture of hard spheres of different diameters d 1 and d 2 . Therefore, if one estimates that a real binary mixture behaves like binary hard-sphere fluid, the variational method consists of minimizing F HS (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 thermodynamic perturbation theory (27). I would like to express my acknowledgement to Jean-Marc Bomont for its stimulating discussions. Thermodynamics is one of the most exciting branches of physical chemistry which has greatly contributed to the modern science. Being concentrated on a wide range of applications of thermodynamics, this book gathers a series of contributions by the finest scientists in the world, gathered in an orderly manner. It can be used in post-graduate courses for students and as a reference book, as it is written in a language pleasing to the reader. It can also serve as a reference material for researchers to whom the thermodynamics is one of the area of interest.