J. Phys. Chem. C 2007, 111, 15957-15961
15957
Thermodynamic Perturbation Theory and Phase Diagram in Simple Fluids† M. Robles and M. Lo´ pez de Haro* Centro de InVestigacio´ n en Energı´a, UniVersidad Nacional Auto´ noma de Me´ xico (U.N.A.M.), Temixco, Morelos 62580, Me´ xico ReceiVed: May 15, 2007; In Final Form: August 20, 2007
Using the lowest-order liquid-state thermodynamic perturbation theory and the hard-sphere fluid as the reference system, the phase diagram of attractive hard-core Yukawa fluids is computed. The results are presented in the reduced temperature versus reduced density, reduced pressure versus reduced temperature, and reduced pressure versus reduced density planes. A comparison with available simulation data and a discussion of the merits and limitations of this approach are also provided.
1. Introduction It is well known1 that two key objectives of the statistical physics of liquids are, on one hand, to understand why and how certain particular phases are stable in well-defined intervals of density and temperature and, on the other hand, to relate the stability, structure, and dynamical properties of such phases with the size and form of the atoms and molecules that constitute the liquid and with the law that governs the intermolecular interaction. In general, these two objectives are very difficult to achieve, and the results obtained so far rely on a variety of simplifications and approximations. Notable among these simplifications and approximations is the assumption of pairwise additive and central forces for the intermolecular interactions. Yet, even for the simplest 3D model that includes this assumption, that of a hard-sphere (HS) fluid, the actual determination of the corresponding phase diagram has defied analytical efforts and is known mainly through computer simulation studies. Nonetheless, for this thoroughly studied system accurate (analytical) approximations exist for the equation of state2 and the structural properties in the fluid region,3,4 and these approximations have served as a starting point for the study of the thermodynamic and structural properties of real fluids. A well-established fact is that in real gases and liquids at high temperature the state and thermodynamic properties are determined almost entirely by the repulsive forces among molecules. At lower temperatures, attractive forces become significant, but even in this case, they affect the configuration of the system at moderate and high densities very little. These findings are taken into account in the application of the perturbation theory of fluids,5-8 where hard-core fluids are used as the reference systems in the computation of the thermodynamic and structural properties of real fluids. However, successful results using perturbation theory are somewhat limited by the fact that, as mentioned above, there are in general no exact (analytical) expressions for the thermodynamic and structural properties of the reference systems that are in principle required in the calculations. The aim of this article is to show that a rather simple perturbation theory treatment with the HS fluid as the reference †
Part of the “Keith E. Gubbins Festschrift”. * Corresponding author. E-mail:
[email protected].
system can yield a reasonably accurate picture for the phase diagram of a simple fluid, in this case, an attractive hard-core Yukawa (HCY) fluid. The choice of this model is due to the fact that it embodies the essential features required for the study of both the liquid state and liquid-gas criticality, namely, a highly repulsive core plus an attractive well, as well as the flexibility of allowing for a variable range of attraction so that, in the appropriate limits, it may serve to model either a LennardJones fluid or a colloidal dispersion. As a matter of fact, one could reasonably wonder at this stage why one should bother to carry out such a study because a few successful theoretical approaches to derive the phase behavior of this system (including the mean spherical approximation to the Ornstein-Zernike equation and even some that use the perturbation theory of liquids with the HS system as a reference) are already available in the literature (e.g., refs 9-33). Two reasons underlie the present effort. On the one hand, we want to assess to what extent the compromise between simplicity and accuracy in the application of a perturbation theory approximation favors such an approach. On the other hand, we also wish to illustrate with a further example the usefulness of the approximate expression for the radial distribution function (RDF) of an HS fluid gHS(r) derived within the so-called rational function approximation (RFA) method.4,34-37 To meet the aforementioned goal, the article is organized as follows. In the next section and in order to make the article self-contained, we provide a rather brief account of the liquidstate perturbation theory and its application to an attractive hardcore Yukawa fluid. Section 3 is devoted to a summary of the RFA method. In section 4, we present the results for the phase diagram of attractive HCY fluids. The article is closed in section 5 with some final comments and concluding remarks. 2. Liquid-State Thermodynamic Perturbation Theory The derivation of the thermodynamic and structural properties of a fluid system becomes a rather formidable problem when one wants to deal with realistic intermolecular interactions. For that reason, since the mid-20th century simplifying attempts to (approximately) solve this problem have been devised, among which the perturbation theories of liquids have played a prominent role.5-8 In this instance, the key idea is to express the actual potential in terms of a reference potential (that of the “unperturbed” system) plus a correction term. This in turn
10.1021/jp073754+ CCC: $37.00 © 2007 American Chemical Society Published on Web 10/25/2007
15958 J. Phys. Chem. C, Vol. 111, No. 43, 2007
Robles and Lo´pez de Haro
implies that the thermodynamic and structural properties of the real system may be expressed in terms of those of the reference system which, of course, should be known. In the case of singlecomponent fluids, a natural choice for the reference system is the hard-sphere fluid that has been thoroughly studied, although as mentioned earlier, even for this simple system the thermodynamic and structural properties are known only approximately. Undoubtedly, the availability of the (analytical) HS radial distribution function (RDF) obtained from the solution to the corresponding Percus-Yevick (PY) equation3 represented a major step toward the successful application of the perturbation theory of liquids to more realistic interparticle potentials. However, the lack of thermodynamic consistency between the virial and compressibility routes to the equation of state present in the PY approximation (as well as in other integral equation theories) is a drawback that may question the results derived from its use within a perturbative treatment. Fortunately for our purposes, another analytical approximation for the RDF of the HS fluid, which avoids the thermodynamic consistency problem, has been more recently derived,4,34-37 and we will describe it in the next section. Suffice it to say here that this development will allow us to examine the phase diagram of the attractive HCY fluid in a rather simple analytical way. Let us now consider a system defined by a pair interaction potential φ(r) split into a known (reference) part φ0(r) and a perturbation part φ1(r). The usual perturbative expansion5 for the (excess) Helmholtz free energy Aex to first order in β ≡ 1/kBT, with T being the absolute temperature and kB being the Boltzmann constant, leads to
Aex Aex 0 ) + 2πFβ NkBT NkBT
∫0∞ φ1(r) g0(r) r2 dr + O(β2)
(1)
where F is the density and Aex 0 and g0(r) are the (excess) free energy and the radial distribution function of the reference system, respectively, and N is the number of particles. If the reference system is taken to be the HS fluid, then of course η Aex HS/NkBT ) ∫0 (ZHS - 1)/η′dη′ and g0(r) ) gHS(r), where ZHS ≡ 1 + 4ηgHS(d+) is the compressibility factor of the HS fluid, η ) (π/6)Fd3 is the packing fraction, d is the HS diameter, and gHS(d+) is the contact value of the HS RDF. Note that given φ(r), in order to use the liquid-state perturbation theory to obtain an explicit approximation to the free energy of the system one has to choose one of the available expressions for ZHS and consider an approximation to gHS(r) for substitution into eq 1. However, depending on the form of the perturbation part of the pair interaction potential, simplified schemes that use the HS fluid as the reference fluid have been introduced in which one usually takes an effective (in general, density- and temperature-dependent) diameter. In ref 40, we have provided a unifying framework for all of the classical perturbative schemes (Barker-Henderson,6 Mansoori-Canfield and Rasaiah-Stell,7 and Weeks-Chandler-Andersen8) considering the approximation to gHS(r) derived with the RFA method (to be briefly described below) and have illustrated it for the case of a Lennard-Jones fluid. In this instance, we examine the attractive HCY fluid whose characteristics will now follow. 2.1. Hard-Core Yukawa Fluid. The attractive hard-core Yukawa fluid is a system whose molecules interact via the pair potential
{
∞,
φHCY(r) )
-0
() [
]
{
r e σ0 z(r σ σ ) φ1(r) ) 0 0 -0 exp , r > σ0 r σ0 0,
() [
]
(3)
By substituting eq 3 into eq 1, one gets the lowest-order approximation to the (excess) free energy of the attractive HCY fluid, namely
Aex Aex HCY HS ) NkBT NkBT 2πFβσ00
[
∫σ∞ exp 0
]
z(r - σ0) gHS(r)r dr (4) σ0
or equivalently
Aex Aex HCY HS ) -12ηβ0ezG(z) NkBT NkBT
(5)
where G(z) ≡ G(t ) z) is the Laplace transform of rgHS(r) (with t being the Laplace transform variable) that will be specified later. For the time being, we should only mention that, once z has been fixed, G(z) will be a function only of the packing fraction η. Moreover, in order to get G(t) using the RFA method, the only input required will be ZHS. Therefore, in this instance one gets an explicit analytical result for the free energy of the attractive HCY fluid within the thermodynamic perturbation theory. Furthermore, in this case one may also readily determine the compressibility factor as
ZHCY )
[ ]
ex ∂ AHCY 1+η ∂η NkT
T,N
) ZHS(η) - 12ηβ0ez
∂ (ηG(z)) (6) ∂η
For the sake of completeness, we finally give the explicit form of the equation of state of the attractive HCY fluid in terms of dimensionless variables P* ) Pσ30/0 (with P being the pressure), F* ) Fσ30, and T* ) kBT/0, namely,
P* ) F*T*ZHS
(π6 F*) - 2πF* e ∂F*∂ (F*G(z)) 2 z
(7)
where it is understood that the packing fraction appearing in both ZHS and G(z) must be expressed in terms of F* (i.e., using F* ≡ (6/π)η). 3. Rational Function Approximation Method for the Hard-Sphere Fluid
r e σ0
z(r - σ0) σ0 exp , r > σ0 r σ0
where σ0 is the hard-core diameter, 0 is the depth of the potential well at r ) σ0, and z is the inverse range parameter. If z = 2, then the properties of this fluid are similar to those of the Lennard-Jones fluid. Clearly, if z f 0 then this potential leads to Coulombic interaction, and in the limit of z f ∞, it reduces to the hard-sphere interaction where the diameter of the spheres is σ0. Also, taking the limits z f ∞ and 0 f ∞ in such a way that the second virial coefficient remains finite, the attractive HCY fluid leads to a version of the adhesive hardsphere model.38,39 We now rewrite φHCY(r) in the form required by the liquidstate perturbation theory as φHCY(r) ) φ0(r) + φ1(r), where the reference potential φ0 will be a hard-sphere pair potential (where the effective diameter d eff of the HS reference system will simply coincide with the hard-core diameter, namely, d eff ≡ σ0) and the perturbed part φ1 is given by
(2)
In this section, we recall the key aspects of the RFA method.4,34-37 The purpose of this analytical-algebraic approach
Thermodynamic Perturbation Theory, Phase Diagram
J. Phys. Chem. C, Vol. 111, No. 43, 2007 15959
is to determine (approximate) reliable expressions for the RDF gHS(r) of a hard-sphere fluid. The main idea of the RFA approach is to retain the form given by Wertheim’s solution to the Percus-Yevick equation3 and use a rational function assumption (hence the name of the procedure) for a function related to the Laplace transform G(t) ≡ L [rgHS(r)] (where L denotes the Laplace transform operator), namely, one takes G(t) to be given by
G(t) )
1 t 12η 1 - etΦ(t)
(8)
Here, Φ(t) is a rational function of the form
1 + S1t + S2t2 + S3t3 + S4t4
Φ(t) )
(9)
1 + L1t + L2t2
S2 ) S3 )
1 η + 12ηL2 + 2 - 24ηS4 2 2η + 1
(10)
L2 ) -3(ZHS - 1)S4
(11)
3 -1 + 4L2 - 8S4 S1 ) η 2 2η + 1
(12)
1 -η + 8ηL2 + 1 - 2L2 - 24ηS4 2 2η + 1
(13)
2 2 1 2η - η + 12ηL2(η - 1) - 1 - 72η S4 (14) 12 (2η + 1)η
S4 )
[ [
(
ZHS - 1/3 χHS 1-η 1 1 + -1 ZHS - ZPY χPY 36η(ZHS - 1/3)
)] ] 1/2
(15)
where ZPY ) (1 + 2η + 3η2)/(1 - η)2 and χPY ) (1 - η)4(1 + 2η)2. Note that the only input required is an expression for ZHS. The radial distribution function, where r is expressed in units of the diameter d, is then given by
gHS(r) )
1
∞
∑ φn(r - n) Θ(r - n)
12ηr n)1
(16)
with Θ(r) being the Heaviside step function and φn(r) ) -L -1[t[Φ(t)]-n]. Explicitly, using the residue theorem, 4
φn(r) ) where
n
e ∑ ∑ i)1 m)1 tir
A(i) mn (n - m)!(m - 1)!
tfti
() [ ] d dt
m-1
t
(t - ti) Φ(t)
n
(18)
with ti being the roots of 1 + S1t + S2t2 + S3t3 + S4t4 ) 0. To guarantee that gHS(r) vanishes inside the core as well as being definitely positive, S4 must be negative. However, according to the form of S4 (cf. eq 15) it may well happen that once ZHS has been chosen there exists a certain packing fraction above which S4 is no longer negative. This was interpreted in refs 35 and 41 as an indication that, for packing fraction η0 where S4 vanishes, the system ceases to be a fluid and a glass transition in the HS fluid occurs. Using this fact, in previous work we have used thermodynamic perturbation theory to study the glass transition both in a Lennard-Jones and in an HCY fluid.42,43 However, this issue will not be discussed further here. 4. Results
Coefficients Si and Li are algebraic functions of η determined by imposing the following two physical restrictions on the RDF: 1. The first integral moment of hHS(r) ) gHS(r) - 1 is welldefined and nonzero. 2. The second integral moment of hHS(r) must guarantee the thermodynamic consistency of the compressibility factor ZHS ) P/(FkBT) and the isothermal susceptibility χHS ) (d(ηZHS)/ dη)-1. After some cumbersome but straightforward algebra (for details, see refs 4 and 34-36), these two conditions lead to
L1 )
A(i) mn ) lim
rn - m
(17)
To proceed and perform actual calculations, one needs a choice for ZHS. Because this quantity is unfortunately not known exactly, as mentioned earlier we will consider the one that follows from the popular Carnahan-Starling equation of state given by2
ZHS(η) )
1 + η + η2 - η3 (1 - η)3
(19)
which is known to be rather accurate in the stable fluid region. For the computation of the gas-liquid coexistence curves, we performed the usual Maxwell construction using eqs 7-15 and 19. The freezing line was determined from the condition F/fr ) (6/π)ηfr where ηfr ) 0.494 corresponds to the packing fraction of the freezing transition in the HS fluid, and the melting line was determined from the condition F/m ) (6/π)ηm where ηm ) 0.545 corresponds to the packing fraction of the melting transition in the HS fluid. Hence, in the F* - T* plane the freezing and melting lines are simply lines parallel to the T* axis whose abscissas are F/fr and F/m, respectively. This result is due to the fact that we are restricting ourselves to the firstorder term in the thermodynamic perturbation expansion. Of course one should expect the finite Yukawa part of the potential to have an influence on the temperature dependence of the freezing and melting transitions, but this is not captured by the present approximation. As far as the location of the critical point is concerned, it was made through the usual route (∂P*/∂F*) ) 0 and (∂2P*/∂F*2) ) 0. Spinodal curves were derived from the maxima and minima of the unstable branch of the reduced pressure versus reduced density isotherms. Finally, the triple point was determined from the intersection of the gas-liquid coexistence curve and the freezing line. Some illustrative results of our computations are displayed graphically in Figures 1-3, where we have chosen to represent the phase behavior of the attractive HCY fluids in the three different thermodynamic planes, namely, the T* - F*, P* T*, and P* - F* planes. In Figure 4, we have represented the critical and triple-point temperatures as functions of the inverse range z-1. Furthermore, in Table 1 we present a comparison of our values of the critical parameters for different z values with those of the theory by Shukla27 and of the self-consistent Ornstein-Zernike (SCOZA) theory.22 As clearly seen in the Figures and in the Table, the qualitative trends observed in the simulations and also obtained with other approaches are indeed reproduced. In particular, for the z dependence of the critical temperature and of the triple-point temperature and for the gas-liquid coexistence in the P* - T*
15960 J. Phys. Chem. C, Vol. 111, No. 43, 2007
Figure 1. Gas-liquid coexistence lines of the attractive hard-core Yukawa fluid in the reduced temperature vs reduced density plane for various values of z as obtained with our perturbation approach using the Carnahan-Starling equation of state for ZHS. Also included are the (vertical) freezing and melting lines and simulation data of Shukla27 (9) and density-functional calculations for z ) 2 (().18 The critical points have been indicated by crosses and have been joined with a dotted line to give insight into the dependence of their location as a function of z. The dashed lines correspond to the spinodal curves derived from the present approach.
Robles and Lo´pez de Haro
Figure 4. Critical temperature and triple-point temperature as functions of z-1. The dotted and continuous lines are the results of the theory by Shukla27 for the critical temperature and the triple-point temperature, respectively. The triangles are simulation data from the literature.19,20,25 The squares are results of the present approach. For comparison, we have also included the critical temperature dependence on z-1 as obtained in the mean spherical approximation via the compressibility route13 (---). Notice the great deviation of this latter both from the results of Shukla27 and from the simulation data, especially for high z. Pentagons (
