3716
J. Phys. Chem. B 2007, 111, 3716-3726
A Perturbation Method for the Ornstein-Zernike Equation and the Generic van der Waals Equation of State for a Square Well Potential Model Byung Chan Eu* and Yuan Qin Department of Chemistry, McGill UniVersity, 801 Sherbrooke St. West, Montreal, Quebec H3A 2K6, Canada ReceiVed: December 15, 2006; In Final Form: February 13, 2007
We calculate the generic van der Waals parameters A and B for a square well model by means of a perturbation theory. To calculate the pair distribution function or the cavity function necessary for the calculation of A and B, we have used the Percus-Yevick integral equation, which is put into an equivalent form by means of the Wiener-Hopf method. This latter method produces a pair of integral equations, which are solved by a perturbation method treating the Mayer function or the well width or the functions in the square well region exterior to the hard core as the perturbation. In the end, the Mayer function times the well width is identified as the perturbation parameter in the present method. In this sense, the present perturbation method is distinct from the existing thermodynamic perturbation theory, which expands the Helmholtz free energy in a perturbation series with the inverse temperature treated as an expansion parameter. The generic van der Waals parameters are explicitly calculated in analytic form as functions of reduced temperature and density. The van der Waals parameters are recovered from them in the limits of vanishing density and high temperature. The equation of state thus obtained is tested against Monte Carlo simulation results and found reliable, provided that the temperature is in the supercritical regime. By scaling the packing fraction with a temperature-dependent hard core, it is suggested to construct an equation of state for fluids with a temperature-dependent hard core that mimicks a soft core repulsive force on the basis of the equation of state derived for the square well model.
I. Introduction The van der Waals equation of state1 describes the thermodynamic properties in a qualitatively correct manner throughout temperature and density ranges although its prediction of the critical properties of fluids is considerably off the experimental values and other thermodynamic properties are not of sufficient accuracy from the current scientific standards. It therefore has spawned numerous ways of modifying its basic ideas in the past2-9 in efforts to understand fluid properties. Such refinements have been mostly phenomenological. On the basis of the virial form of equation of state, it was shown in a series of papers10-12 that the van der Waals equation of state could be generalized to a form formally resembling it without an approximation. We called it the generic van der Waals (GvdW) equation of state, which indicates it is a canonical form of equation of state. The parameters therein are called the generic van der Waals (GvdW) parameters. They are generally functions of temperature and density that tend to the van der Waals constants, as the density diminishes to zero and the temperature increases to infinity. The GvdW equation of state has been applied to study phenomenologically the critical isotherm12 and to successfully calculate various transport properties13,14 of liquids from the molecular theory standpoint by means of Monte Carlo simulations or a numerical solution method15-18 for the Percus-Yevick integral equation for the pair correlation function. The GvdW equation of state for a square well potential10 in particular is simply given in terms of the cavity functions at the contact diameter of the hard sphere and at the outer position of the well. Thus, the form of the GvdW equation of state for a square well potential is rather suggestive and tempting, but the cavity functions at the points mentioned are not known analytically in closed form. They, however, can be calculated
numerically as a function of density and temperature by a method of simulations, but there are still practical needs to have the GvdW parameters in closed forms as functions of temperature and density for real fluids. Therefore, some sort of analytic form is much desired for them. To meet this kind of necessity, we have developed empirical models for GvdW parameters A and B, especially for the subcritical regime,11 and studied the critical isotherm12 as a function of density, which successfully correlates with experiment. In this article, we calculate the GvdW parameters for supercritcal fluids by using a square well potential model by means of a perturbation method for the Percus-Yevick (PY) integral equation. Admittedly, square well potentials are not fully realistic for a thermodynamic description of real spherical fluids, but we find that with the model it is possible to coax analytical forms for A and B to emerge in terms of density and temperature, if we use the PY integral equation and solve the equations by a perturbation method treating either the well width or the Mayer function [exp(/kBT) - 1] ( is the well depth, kB is the Boltzmann constant, and T is the temperature) or functions in the square well region exterior to the hard core as the perturbation parameter. In fact, the perturbation terms invariably come under integral operators in the integral equations. Therefore, the present perturbation method may be said to be a variant of the Neumann series method for Volterra-type integral equations,19 which PY integral equations are. We then calculate the GvdW parameters B and A by the perturbation solutions. Even at their lowest order, we find A and B are already considerably more generalized than the original van der Waals parameters and those existing in the literature.6,7 In 1964, Longuet-Higgins and Widom6 proposed what amounted to an empirical modification of the van der Waals
10.1021/jp068641p CCC: $37.00 © 2007 American Chemical Society Published on Web 03/20/2007
Square Well Potential Model equation of state in which the repulsive force contribution to the equation of state was identified with the Alder and Wainwright simulation result20 whereas the attractive force contribution was taken the same as the original van der Waals parameter a. With such a modification, thermodynamic properties except for certain second derivatives of the free energy near the triple point were found to agree well with experiment. Such an idea was taken up by Guggenheim,7 who showed that one could construct similar equations of state with the repulsive contribution identified with the scaled particle theory result21 or the hard sphere result22-24 obtained from the PY integral equation. More recently, the same idea was taken up for studies of inherent structures25,26 in connection with glass transition and void distribution27 by a computer (Monte Carlo) simulation method. In these empirical forms of equation of state the van der Waals parameter a is also assumed to be a pure constant independent of density and temperature, as was for the original van der Waals theory. What we obtain in this work is a statistical mechanical derivation of analytic forms for the equation of state representing those mentioned, albeit by using a perturbation method and thus being approximate. In the literature, thermodynamic properties have been studied by perturbation theory in statistical mechanics in the line of approach by Barker and Henderson,28 Weeks, Chandler, and Andersen,29 and numerous other authors after them. However, in the aforementioned class of perturbation theories the Helmholtz free energy is expanded around a suitably chosen reference state Helmholtz free energy, usually that of the hard sphere fluid, which is then calculated by using the pair correlation function for the hard sphere fluid obtained numerically. In the present perturbation theory, although the reference state is still that of the hard sphere fluid, the perturbation solutions are sought for the pair distribution functions after the integral equations for the distribution functions have been suitably recast into coupled integral equations by applying the Wiener-Hopf method for integral equations.14,24,30-34 The theory formulated is capable of producing analytical equations of state and other thermodynamic properties, and the equations of state obtained are sufficiently accurate provided that the temperature is above and away from the critical temperature, as will be shown later. This limitation seems to be mainly rooted in the fact that the integral equation used, namely, the PY integral equation, is not good for description of the critical regime of fluids. Since the integral equations used are amenable to further improvements in that respect and other aspects and also further generalizations to complex fluids, we believe that the present method appears to present a considerable potential for more refined theories of wider applicability and, in that sense, is quite suggestive of a more extensive investigation in the future. At this point, it is also appropriate to mention some existing equations of state developed for associating fluids35 and hard chain molecules.36 The former is based on the statistical mechanical theory originally developed by Wertheim37 for cluster (dimer, trimer, etc.)-forming fluids, namely, dimers, trimers, etc., of hard spheres, whereas the latter generalizes the basic idea of the Flory-Huggins theory38,39 and can give rise to analytic forms for the equation of state for athermal hard chain fluids. Although interesting and suggestive from the standpoint of the present perturbation theory they, however, do not have a direct relevance to the present theory of a simple fluid obeying a square well potential. This paper is organized as follows. In section II, the generic van der Waals equation of state is briefly reviewed. Then the statistical mechanical forms for the GvdW parameters are
J. Phys. Chem. B, Vol. 111, No. 14, 2007 3717 presented for a square well potential model. In section III, the PY integral equation is put into an equivalent pair of integral equations by using the Wiener-Hopf method14,24,30-34 for the square well potential. The pair of integral equations is then solved by a perturbation method.28,29,40 In section IV, with the solutions for the cavity functions we calculate the GvdW parameters as functions of density and temperature to the leading order, that is, B to the zeroth order and A to the first order in the perturbation parameter. Therewith the equations of state are calculated to the first order. In section V, the equations of state are tested in comparison with Monte Carlo (MC) simulation results and found reliable provided that the temperature is in the supercritical regime. Section VI is for discussion and concluding remarks. II. Brief Review of the Generic van der Waals Equation of State A. Generic van der Waals Equation of State. Statistical mechanics yields the virial equation of state in the form34,41
p)
F 2π 2 - F β 3
∫0∞dr r3u′(r) g(r,F,β)
(1)
where β ) 1/kBT, u′(r) ) du/dr, p is the pressure, F is the number density, and g(r,F,β) is the pair correlation function. If the interaction potential u(r) can be separated into the following form
u(r) ) ur(r) for r e σ ) ua(r) for r > σ
(2)
where ur(r) and ua(r) are respectively the repulsive and attractive part of the potentialsthe LJ potential is clearly an example for this kind of potential energysthen the virial equation of state (1) can be cast into the form10
(p + AF2)(1 - BF) ) Fβ-1
(3)
where
A(F,β) )
∫σ∞dr r3u′(r) g(r,F,β)
2π 3
∫
2πβ σ dr r3u′(r) g(r,F,β) 3 0 B(F,β) ) 2πβ σ F 0 dr r3u′(r) g(r,F,β) 13 -
∫
(4)
(5)
Equation 3 is called the GvdW equation of state because the original van der Waals equation of state is a special case of it, and it is a generic and canonical form of equation of state. The parameters A and B are called the GvdW parameters. They are generally dependent on F and T. The point of splitting the virial integral in eq 1, however, is not unique, and recently, it was found42 that there exists a value rm such that
r‡ ) rm/σ < 1 which is the value of r at which the integrand
I ) -r3u′(r) g(r,F,β) is a maximum. If this point is used, then the GvdW parameters are defined by
3718 J. Phys. Chem. B, Vol. 111, No. 14, 2007
A(F,β) )
∫r∞σdr r3u′(r) g(r,F,β)
2π 3
‡
∫
2π r‡σ dr r3u′(r) g(r,F,β) 3 0 B(F,β) ) 2πβ r‡σ 1F 0 dr r3u′(r) g(r,F,β) 3 -
∫
Eu and Qin
(6)
A(F,β) ) a′
(8)
lim
B(F,β) ) b′
(9)
where a′ and b′ are constants independent of F and T. These constants indeed can be identified with the van der Waals parameters as will become evident from the perturbation theory results given later. The GvdW equation of state in these limits is then identified with the van der Waals equation of state1,43
(p + a′F2)(1 - b′F) ) FkBT
(10)
This is the reason for the name given for eq 3. The formulas for A(F,β) and B(F,β) given in eqs 6 and 7 are the statistical mechanical representations for the GvdW parameters, whose limiting forms give rise to the van der Waals parameters. This limiting behavior (8) and (9) is the reason why the van der Waals equation of state cannot be derived exactly from statistical mechanics unless some sort of approximation is made. The GvdW parameters are not possible to obtain in closed forms for realistic potential models since for such potentials the pair correlation functions cannot be computed in closed forms. The hard sphere22-24 and adhesive hard sphere model44 are the only potential models that can afford us an analytic evaluation of the GvdW parameters, but they are not sufficiently realistic. The parameter r‡ was found to be the diameter of a levitation cavity45 in the fluid. For the case of the Lennard-Jones potential, it was found possible to fit it to the mathematical formula42
r (F,T*) ) d0[1 + d1xT*] ‡
-1/6
) - if σ < r < d ) 0 if r > d
lim
Ff0,Tf∞
u(r) ) ∞ if r < σ
(7)
By using the density and temperature dependence of g(r) empirically or numerically or approximately known in the literature, it is possible to deduce the limiting behavior10 of A(F,β) and B(F,β): Ff0,Tf∞
curtailed and thus of rather short range. Therefore, a square well potential model is not unrealistic in such a situation. The square well potential is mathematically represented by
(11)
where
(14)
where (d - σ) is the width of the wellsthe range of the attractive well. Since it is convenient to work with reduced dimensionless variables, we define the reduced variables appearing in the theory presented below:
π π x ) r/σ; γ ) d/σ; η ) σ3F; F* ) Fσ3, V0 ) σ3 6 6 β* ) β; * ) /kBT ) 1/T*; p* ) pV0/ A*(η, T*) ) A(F,T)/V0; B*(η,T*) ) B(F, T)/V0; K* ) K/σ2 where K(r) is the Wiener-Hopf factorization function that will appear presently in connection with the PY integral equation. Its precise meaning will become clear when we deal with the integral equation according to the Wiener-Hopf method.14,24,30-34 With this system of reduced variables, the GvdW parameters for the square well potential model are given by the formulas
A*(η,T*) )
4(e* - 1) 3 [γ y(γ,η,T*) - y(1,η,T*)] (15) *
B*(η,T*) )
4y(1,η,T) 1 + 4ηy(1,η,T)
(16)
In these expressions, y(γ,η,T*) and y(1,η,T*) are the cavity functions defined by
y(x) ) exp[βu(x)]g(x)
(17)
which are evaluated respectively at x ) γ and x ) 1. They should be obtained by solving the integral equation for the pair correlation function g(x), for example, the PY integral equation. III. The Integral Equations and Their Solutions
d0 ) 1.10956, d1 ) 0.76531 + 0.07737F* 0.009937F*2 + 0.10660F*3 (0.2 j F* ) Fσ3 j 1.1) (12) The temperature dependence of this parameter originates from the fact that the repulsive potential is soft, so that the effective contact diameter of a pair of particles gets smaller as the temperature increases. Since the density dependence of r‡(F,T*) is rather weak, we will find it convenient to take -1/6
r (T*) ) [1 + 0.76531xT*] ‡
A. The PY Integral Equations. There are various integral equations47,48 that can be used to determine the pair correlation function. The PY integral equation is one example, which is most extensively studied and used in the literature. It is obtained from the Ornstein-Zernike (OZ) integral equation41
∫
h(x) ) c(x) + F* dx′ c(|x1 - x′|)h(|x′- x2|) (F* ) Fσ3, x ) |x1 - x2|) (18)
(13)
as a reasonable approximation (d0 ≡ 1), taking densityindependent parameters. It should be noted that a form similar to the one given in eq 13 was originally used by Boltzmann.46 We may make use of this variable contact diameter to construct an empirical equation of state as suggested in section VI. B. Square Well Potential Model. In high-density liquids in which particles are closely packed, the repulsive interaction is almost of hard spheres and the attractive interaction is much
in the reduced variables used. Here h(x) is the total correlation function h(x) ) g(x) - 1 and c(x) is the direct correlation function. If the OZ equation is closed by the PY closure
c(x) ) f(x) y(x)
(19)
f(x) ) exp[-βu(x)] - 1
(20)
where
Square Well Potential Model
J. Phys. Chem. B, Vol. 111, No. 14, 2007 3719
is the Mayer function, then the PY integral equation arises from the OZ equation. It can be solved by a variety of methods. In the Wiener-Hopf method of solution for the integral equations, the PY integral equation is Fourier transformed and the Fourier transform C ˆ (k) of the direct correlation function c(x)
h(x) ) -1 for x < 1
C ˆ (k) )
∫0 dx C(x) cos kx ∞
C(x) ) 24η
∫x∞ds sc(s)
) e*y(x) - 1 for 1 < x < γ ) y(x) - 1 for x g γ
In the case of the square well potential, it is convenient to denote K*(x) in different intervals of x by different symbols. Thus, we define
K/1(x) ) K*(x) for 0 < x < 1 ) 0 for x g 1
is factorized to a product of real functions in the form
(21)
K/2(x)
The factorization function K ˆ (k) is the Fourier transform of the real function K*(x) defined in x space through the relation
K/3(x)
1-C ˆ (k) ) K ˆ (k) K ˆ (-k)
1-K ˆ (k) ) 12η
∫-∞∞dx K*(x) exp(ikx)
(22)
It should satisfy certain conditions such as regularity, the absence of zeros on the real axis, and a vanishing limit of |1 - K ˆ (k)|, namely, |1 - K ˆ (k)| ∼ O(1/k), as |k| f ∞. The inverse relation of eq 22 is
12ηK*(x) )
1 2π
∫-∞∞dk [1 - Kˆ (k)] exp(-ikx)
(23)
It then is possible to show that the PY integral equation can be written as a pair of integral equations for the direct and total correlation functions c(x) and h(x), respectively, given in terms of K*(x):
xc(x) ) -K*′(x) + 12η xh(x) ) -K*′(x) + 12η
∫x∞ds K*′(s) K*(s - x)
∫0∞ds(x - s) h(|x - s|) K*(s)
(24) (25)
where the prime denotes differentiation with respect to x. For a finite ranged potential the upper limit of the convolution integral of the set is finite, say, γ for the case of the square well potential mentioned. Derivation of the set of eqs 24 and 25 is well documented, and the reader is referred to the literature.24,34 Since there are three unknown functions in the equations presented, another equation is required to solve them. Such an equation is provided by a closure relation. With the closure relation (19), eqs 24 and 25 are a closed pair of integral equations for the cavity function y(x) and the factorization function K*(x). It should be noted that
K*(x) ) 0 for x < 0
) K*(x) ) 0 for x g γ
(28)
such that by continuity of K*(x)
K/1(1) ) K/2(1), K/2(γ) ) K/3(γ)
(29)
Henceforth, the asterisk will be omitted from K/1(x) and K/2(x) for brevity of notation. The factorization functions appearing below therefore should be understood for dimensionless functions. We also distinguish the cavity function in different intervals:
y1(x) ) y(x) for 0 < x < 1 y2(x) ) y(x) for 1 < x < γ y3(x) ) y(x) for x g γ
(30)
which are also continuous at x ) 1 and x ) γ. Now, there are two different sets of integral equations for y1(x), y2(x), K1(x), and K2(x) in the range 0 < x < γ:
∫x1ds K′1(s) K1(s - x) + 12ηOy1(K1,K2|x)
xy1(x) ) K′1 - 12η
(31) K′1(x) ) x - 12η
∫01ds(x - s) K1(s) + 12ηOK1(K2,y2|x)
xy2(x) ) x - 12η
(32)
1 ds(x - s) K1(s) + 12ηOy2(K1,K2,y2|x) ∫x-1
(33) 1 ds(x - s) K1(s) + ∫x-1
12ηOK2(K1,K2,y2|x) (34)
if the potential energy is finite ranged and vanishes for x > γ. The factorization function K*(x) is a continuous function of x and, in the case of the square well potential, is nonvanishing in the range 0 < x < γ. For the square well potential model subject to the PY closure47
c(x) ) -y(x) for x < 1 ) (e - 1)y(x) for 1 < x < γ *
whereas
) K*(x) for 1 < x < γ ) 0 for x e 1
K′2(x) ) -fx + 12ηf
K*(x) ) 0 for x g γ
) 0 for x g γ
(27)
(26)
where the integral operators consisting of convolution integral operators are defined by
Oy1(K1,K2|x) ) -
∫x1+xds K′2(s) K1(s - x) γ ds K′2(s) K2(s - x) ∫1+x
(35)
∫1γds(x - s) K2(s) γ (1 + f)∫1 ds(s - x) y2(s - x) K2(s)
(36)
OK1(K2, y2|x) ) -
3720 J. Phys. Chem. B, Vol. 111, No. 14, 2007
Eu and Qin
∫0x-1ds(x - s)[y2(x - s) - 1]K1(s) + x-1 γ f∫0 ds(x - s) y2(x - s) K1(s) - ∫1 ds(x - s) K2(s) γ ds K′2(s) K2(s - x) + ∫x1+xds K′2(s) K1(s - x) - ∫1+x γ (1 + f)∫1+xds(x - s) y2(s - x) K2(s) (37)
Oy2(K1, K2, y2|x) )
OK2(K1,K2,y2|x) ) -fOy2(K1,K2,y2|x) + γ ds K′2(s) K2(s - x) ∫x1+xds K′2(s) K1(s - x) + ∫1+x
(38)
Here f ) exp(*) - 1. The integral equation for K3(x) ≡ K(x) in the range x > γ does not exist and that for y3(x) can be obtained easily from eq 25, but since it is not relevant to the question of interest in this work, we shall omit it in this paper. The sets of integral equations presented for the square well cannot be solved in closed form. However, since the hard sphere part can be solved exactly if the attractive part is absent, it suggests a perturbation method based on the hard sphere solution as the reference solution. This is the method we employ in this work. The idea of treating the hard sphere part as the zeroth order solution is similar to that of the thermodynamic perturbation theories28,29 in which the Helmholtz free energy is expanded49 in a perturbation series by treating the attractive contribution as a perturbation. The convolution integral operators (35)-(38) invariably involve functions in region II (1 < x < γ) and thus pertain to the square well (i.e., the attractive potential). On inspecting eqs 31-34, we observe that, except for the K′2(x) term the terms in region II are either multiplied by f associated with the attractive potential or integrated over region II (namely, integral operators) whereas the terms in region I (0 < x < 1) are not. Thus, it is reasonable to treat the integral operator terms in region II as perturbations to the hard sphere contributions arising from region I. On the basis of this observation, we will multiply a perturbation parameter λsthe book-keeping parameter for the moments to the integral operators and develop a perturbation theory of solution for the integral equations presented. The perturbation parameter λ will be eventually set equal to unity in the end, but the end results suggest that it may be identified with f(γ - 1) ≡ fδ since the well width δ appear invariably multiplied by the factor f. This is not surprising because the attractive potential is in essence treated as the perturbation. On multiplying the parameter λ to the integral operators, we obtain the governing integral equations in the forms
xy1(x) ) K′1 - 12η
1
12ληOy1(K1,K2|x) (39)
∫01ds(x - s) K1(s) + 12ληOK1(K2,y2|x) (40)
xy2(x) ) x - 12η
(1) 2 (2) yi(x) ) y(0) i (x) + λyi (x) + λ yi (x) + ... (i ) 1, 2)
1 ds(x - s) K1(s) + ∫x-1
These expansions are physically sensible if λ is smaller than 1 and thus expected to yield reasonable results for large hard core particles with a short-ranged attractive potential. This situation is expected to be generally realizable for liquids because the attractive potential is effectively curtailed owing to the proximity of repulsive particles in liquid states. Inserting these expansions into eqs 39, 40, 41, and 42 and equating the terms of like power of λ, we obtain hierarchies of (j) equations for the expansion coefficients y(j) i and Ki (x):
K(0)′ 1 (x) ) x - 12η
∫01ds(x - s) K(0)1 (s)
1 ds(x - s) K1(s) + ∫x-1
12ληOK2(K1,K2,y2|x) (42)
It is important to keep in mind that eqs 39 and 40 are for x < 1 whereas eqs 41 and 42 are for 1 < x < γ. These coupled
(45)
∫01ds(x - s) K(1)1 (s) + 12ηOK1(K(0)2 ,y(0)2 |x)
K(1)′ 1 (x) ) -12η
(46) l, K(0)′ 2 (x) ) 0 K(1)′ 2 (x) ) -fx + 12ηf
(47)
1 ds(x - s) K(0) ∫x-1 1 (s) + (0) (0) 12ηOK2(K(0) 1 ,K2 , y2 |x) (48)
l, (0)′ xy(0) 1 (x) ) K1 - 12η
(0) ∫x1ds K(0)′ 1 (s) K1 (s - x)
(49)
(1) ∫x1ds[K(0)′ 1 (s) K1 (s - x) +
(0) (0) (0) K(1)′ 1 (s) K1 (s - x)] + 12ηOy1(K1 ,K2 |x) (50)
l, xy(0) 2 (x) ) x - 12η xy(1) 2 (x) ) x - 12η
1 ds(x - s) K(0) ∫x-1 1 (s)
(51)
1 ds(x - s) K(1) ∫x-1 1 (s) +
12ληOy2(K1,K2,y2|x) (41)
K′2(x) ) -λfx + 12ληf
(43)
(1) 2 (2) Ki(x) ) K(0) i (x) + λKi (x) + λ Ki (x) + ... (i ) 1, 2) (44)
(1)′ xy(1) 1 (x) ) K1 - 12η
∫x ds K′1(s) K1(s - x) +
K′1(x) ) x - 12η
integral eqs 39, 40, 41, and 42 will be solved to obtain the cavity functions at x ) 1 and x ) γ. Before proceeding further, we note that these integral equations are of Volterra type,19 which may be solved in a Neumann series. A perturbation theory for such integral equations is in fact a method of generating Neumann series as solutions for them, in which the inhomogeneous terms are taken with the hard sphere solution, especially because the dynamics of particles in liquids is dominated by that of hard spheres and the hard sphere problem can be exactly solved in the case of the PY integral equation. We are exploiting this latter feature of the problem to obtain analytic results for the GvdW parameters and the equation of state for the square well model to the order of solutions considered. B. Perturbation Solutions. We expand the solutions as follows:
(0) (0) 12ηOy2(K(0) 1 ,K2 ,y2 |x) (52)
l. The hierarchies of the perturbation equations presented, (45)(52), can be solved without difficulty order by order. As a matter of fact, to calculate B* to O(1) and A* to O(fδ), it is sufficient (0) (0) to obtain only K(0) 1 , y1 , and y2 . In this regard, it should be
Square Well Potential Model
J. Phys. Chem. B, Vol. 111, No. 14, 2007 3721
remarked that in the thermodynamic perturbation theory the free energy is calculated to O(fδ) or the equivalent. For the aforementioned purpose we only need to solve eqs 45, 49, and 51. Equation 45 is not coupled to the rest of the hierarchies and yields the well-known hard sphere solution14,24 for the factorization function. This solution is used to calculate the convolution integral in eq 49 and thus y(0) 1 , which also gives the hard sphere cavity function. Using the continuity condition
Note that we now have set λ ) 1, and the product fδ may now be regarded as the perturbation parameter. With A* and B* in eqs 56 and 57 the following GvdW equation of state via the virial route is obtained, if we neglect the terms of linear in δ or higher in the square bracket term in eq 56:
(0) K(0) 1 (1) ) K2 (1) ) 0
It should be noted that the second term on the right, namely, the attractive potential contribution, is temperature-dependent, and hence the equation of state is no longer athermal in contrast to the hard sphere equation of state. B. Compressibility Equation of State. On the other hand, since the equation of state via the compressibility route for the square well potential model may be written as
the desired solution for K(0) 1 (x) is easily obtained:
K(0) 1 (x) )
(1 + 2η) 2(1 - η)2
(x2 - 1) -
3η (x - 1) (53) 2(1 - η)2
The procedure to obtain it from eq 45 is well documented in the literature14,24,34 on hard sphere fluids, to which the reader is referred for the details. Using the result for K(0) 1 (x), we now calculate the cavity functions:
y(0) 1 (x) ) xy(0) 2 (x) )
η(2η + 1)2 2(1 - η)4 η(2η + 1) 2(1 - η)
2
x3 -
3η(η + 2)2 2(1 - η)4
x+
(2η + 1)2 (1 - η)4
(54)
∫0ηdη′ η′∫01dx x2y1(x) 24f η ∫ dη′ η′∫1γdx x2y2(x) η 0
the GvdW parameters via the compressibility route are given by
3η(η + 2) 2 3η2 x3 x + 2 (1 - η) (η - 1)2 (6η + 17η2 + 1) 12η2 x(55) 2 (η - 1) (1 - η)2
IV. The GVDW Parameters and Equations of State A. Virial Equation of State. With the solutions obtained (0) for y(0) 1 (x), and y2 (x) it is easy to calculate the GvdW parameters as functions of temperature and density by the virial route: we obtain B* to O(1) and A* to O(fδ). They are as follows:
[
(1 - η + 2η2) (1 - η)2
1+
(4η2 - 7η + 2) 2(2η2 - η + 1)
δ
]
) 12fδT*[1 + η + 3η2 + 5η3 + 7η4 + O(η5)] (56) B*(η,T*) )
2(2 + η) (1 + 2η + 3η2)
where
δ)γ-1
(57)
(58)
p* 24 )1+ ηT* η
x4 -
These lowest order cavity functions can be made use of to calculate the leading order GvdW parameters and thus the equation of state. The higher order cavity functions can be straightforwardly obtained after lengthy calculations by solving eqs 46, 48, 50, and 52. They give rise to rather complicated GvdW parameters and the corresponding equation of state, which does not appear to be practically easy to use. Furthermore, since one of motivations for the present work is to extract relatively simple analytic forms for the equation of state that can be made semiempirical, perhaps, by adding one or two empirical parameters, too complicated forms for the equation of state are not desirable. For this reason we do not present them here.
A*(η,T*) ) 12fδT*
η(1 - η + 2η2) p* 1 + 2η + 3η2 12fδ ) ηT* (1 - η)2 (1 - η)2
A* )
B* )
∫0ηdη′ η′∫1γdx x2y2(x,η′)
24fT* η2 24 η2 1+
(59)
∫0ηdη′ η′∫01dx x2y1(x, η′) 24 η
(60)
∫0ηdη′ η′∫01dx x2y1(x,η′)
Using the cavity functions obtained, we find the GvdW parameters by the compressibility route in the forms
A* )
[
12fδT* (4 - η) 4 + ln (1 - η) η 1-η η
]
5 11 7 ) 12fδT* 1 + η + 2η2 + η3 + η4 + O(η6) 3 5 3
[
]
B* )
(4 - 2η + η2)η
(61)
(62)
(1 + η + η2)
These are different from those by the virial route for the same reason as for the difference in the equations of state by the two routes. On substitution of these formulas into the GvdW equation of state (3), the compressibility equation of state to O(δ) is obtained:
[
]
12(4 - η) 48 1 + η + η2 p* - fδ ) + ln(1 - η) 3 ηT* 1-η η (1 - η)
(63)
It is known in the case of hard sphere fluids that a linear combination of the virial and compressibility equations of state yields the Carnahan-Starling equation of state:34,50
p* 1 p* 2 p* 1 + η + η2 - η3 ) + ) (64) ηT* 3 ηT* virial 3 ηT* comp (1 - η)3
( )
( )
Therefore, it is tempting to use the same linear combination to obtain the equation of state for the hard sphere part with the
3722 J. Phys. Chem. B, Vol. 111, No. 14, 2007
Eu and Qin
Carnahan-Starling form. We thereby obtain for the virial equation of state to O(fδ)
(1 - η + 2η2) p* 1 + η + η 2 - η3 ) 12fδη (65) ηT* (1 - η)3 (1 - η)2 whereas for the compressibility equation of state
p* 1 + η + η 2 - η3 ) ηT* (1 - η)3
[
fδ )
]
12(4 - η) 48 + ln(1 - η) (66) 1-η η
1 + η + η 2 - η3 (1 - η)3 5 11 12fδη 1 + η + 2η2 + η3 + ... (67) 3 5
(
)
These in fact will be used for numerical comparison and examining the empirical GvdW parameters employed in ref 12. V. Assessment of the Theoretical Results A. Numerical Comparison with Simulation Results. By treating the equations of state (65) and (66) as semiempirical, it may be possible to apply them to experimental systems to study their thermodynamic properties. However, since such application is not the immediate aim of this work, it is not attempted here. Instead, to ascertain their utility, we test them in comparison with MC simulation results for the square well potential. The numerical test in the following will be performed with the equations of state (65) and (66). Therefore, the difference between the two forms of equation of state is in the attractive contributions only. They have been found to slightly improve the accuracy over the original perturbation theory forms, eq 58 and eq 63, respectively. B. Description of the Monte Carlo Simulation Method. Here the MC simulation method51,52 used for generating the numerical equations of state is briefly described. 256 particles, subjected to the square well potential of interaction, are placed in a cubic box and simulated. A face-centered cubic lattice is used as the starting configuration, and the Metropolis algorithm is used to generate the NVT ensemble distribution of the particles. The maximum step size is adjusted, so that 50% of the attempted moves are accepted. The first 2 × 105 MC cycles are rejected to eliminate initial configuration memory, and the subsequent 1 × 106 cycles are collected to calculate the radial distribution functions. The pressure is calculated from the virial equation of state
p*β* ) 1 + 4ηy(σ) + 4fη[y(σ) - γ3y(γ)] η
(68)
obtained from eq 1 for the square well. The pair correlation functions y(σ) and y(γ) for this equation are obtained by extrapolating or interpolating the radial distribution functions. C. Numerical Comparison of Theoretical and Simulation Results. In the present theory there appear two distinct parameters, f ) exp(1/T*) - 1 and δ ) γ - 1. The former may be, of course, regarded as a measure of the potential well depth, and the latter is the range of the well. In the present theory, since not only the radius of convergence of the perturbation expansion is not known but also the complicated higher order results are not desirable for practical applications,
Figure 1. Comparison of theoretical results with MC simulations. δ ) 0.05. The symbols and the numerals for solid and broken curves are MC simulation values and the perturbation theory results, respectively: 1 and open circles (O) for T* ) 4.0; 2 and triangle (4) for T* ) 1.0; and 3 and diamond (]) for T* ) 0.6. The solid and broken curves are the perturbation theory results obtained from the virial and compressibility equation of state, respectively.
we test eqs 65 and 66 in comparison with MC simulation results for various values of δ as well as different reduced temperatures that are believably within the supercritical regime. (The equation of state (65), for example, does not yield the critical point in the fluid regimesnamely, it yields unphysical critical points but if the attractive part is expanded, then there exists a critical point in the fluid regime, although quite off the experimental value for a simple fluid, as shown later.) In Figure 1, the reduced equations of state are plotted against the packing fraction η ) πσ3F/6 for three different reduced temperatures: T* ) 4.0, 1.0, and 0.6 in the case of δ ) 0.05. The solid and broken curves distinguished by numerals are respectively the virial equation of state and the compressibility equation of state. The symbols are the MC simulation results. The meanings of the numerals and symbols are as follows: 1 and O for T* ) 4.0, 2 and 4 for T* ) 1.0, and 3 and ] for T* ) 0.6. In this case of δ ) 0.05, the model is, in fact, almost a sticky hard sphere, but there appear considerable deviations between the simulation and theoretical results, although the theory predicts a qualitatively correct behavior with respect to η. As will be shown, comparison gets improved in quality as δ gets larger in value. There are two reasons that we believe can be attributed to this variance. One that is probably of minor significance is that the MC simulation becomes increasingly difficult as the δ value decreases below a certain point, and it seems reasonable to attribute the deviations observed between theory and simulation values, to some extent, to the difficulty of MC simulations mentioned, which may have larger statistical errors when δ values are too small. The other more significant and major reason is that the temperature (T* ) 0.6) in question is near the critical point of the fluid and the integral equation, namely, the PY integral equation, employed is known to be incapable of properly describing the fluid behavior near and at the critical point. We believe that to describe the critical behavior of fluids the PY integral equation must be significantly improved, but it is not clear at present how that can be achieved for the integral equations to be employed in the present perturbation theory, although there are some clues to the ways to improve it. The investigation into this question should be deferred to the future work. Despite this difficulty in the lower end of the temperature range, it is interesting to note that unlike
Square Well Potential Model
J. Phys. Chem. B, Vol. 111, No. 14, 2007 3723
Figure 2. Comparison of theoretical results with MC simulations. δ ) 0.1. The meanings of the symbols, the numerals, and curves are similar to Figure 1: 1 and open circles (O) for T* ) 4.0; 2 and triangle (4) for T* ) 1.2; and 3 and diamond (]) for T* ) 0.8.
Figure 3. Comparison of theoretical results with MC simulations. δ ) 0.3. The meanings of symbols, the numerals, and curves are similar to Figure 1: 1 and open circles (O) for T* ) 6.0; 2 and triangle (4) for T* ) 2.5; and 3 and diamond (]) for T* ) 1.5.
the case of hard sphere fluids the virial equation of state is invariably of better accuracy than the compressibility equation of state when the attractive contribution is taken into consideration. In Figure 2, a similar comparison is made in the case of δ ) 0.1 and T* ) 4.0, 1.2, 0.8. The meanings of the symbols and numerals are similar to those for Figure 1. For this case of δ, agreement between the simulation result and the theoretical prediction is improved. Again, the virial equation of state is found more accurate than the compressibility equation of state when the attractive contribution is taken into consideration, but in the case of T* ) 0.8 there still is a considerable deviation from the MC simulation result as η increases, indicating that the perturbation theory based on the PY integral equation is getting unreliable as the temperature gets lowered toward the critical point. In Figure 3, similar comparison is continued for the case of δ ) 0.3 and T* ) 6.0, 2.5, 1.5, which shows that for the cases of higher reduced temperature agreement is excellent, and it seems to support the difficulty of simulations in the case of very small δ values mentioned earlier and associated lower degrees of agreement. Nevertheless, at T* ) 1.5 the theory also
does not agree well with the simulation results. We believe that this variance is attributable to the poor quality of the PY integral equation near the critical point which appears near T* ) 0.82 for the value of δ, and it appears that T* ) 1.5 is still too close to T* ) 0.82. In Figure 4, we take δ ) 0.5 and T* ) 10.0, 4.0, and 2.5 and make comparison with simulation data. The compressibility factors are found to be satisfactory as for the previous cases, and at the lowest value of T*, the situation is much improved although not as good as for other temperatures. In Table 1, the numerical values are given for the MC simulation results presented in Figures 1-4. In this tables the second row of the heading is for reduced temperature and the first column is for packing fraction. D. Comparison of Empirical GvdW Parameters and the Perturbation Theory Results. The present work has been largely motivated by the desire to obtain useful analytic forms for the GvdW parameters and understand their empirical formss for example, the quadratic model12ssuccessfully employed for study of critical isotherms of simple fluids. The empirical forms had to be singular functions of density and temperature, so that they give rise to discontinuity in the liquid-vapor coexistence
TABLE 1: MC Simulated Compressibility Factor for Square Well Fluids δ ) 0.05
δ ) 0.1
δ ) 0.2
η
T* ) 4.0
T* ) 1.0
T* ) 0.6
T* ) 4.0
T* ) 1.2
T* ) 0.8
T* ) 4.0
T* ) 2.0
T* ) 1.2
0.10 0.15 0.20 0.25 0.30 0.35 0.40 0.45
1.493 1.852 2.314 2.937 3.759 4.858 6.396 8.607
1.362 1.619 1.931 2.341 2.868 3.601 4.609 6.094
1.173 1.299 1.433 1.640 1.912 2.294 2.860 3.647
1.465 1.923 2.231 2.805 3.564 4.621 6.061 8.183
1.287 1.481 1.738 2.089 2.535 3.169 4.101 5.509
1.110 1.187 1.326 1.498 1.734 2.135 2.711 3.643
1.402 1.699 2.078 2.596 3.294 4.312 5.740 7.868
1.275 1.468 1.734 2.114 2.636 3.413 4.591 6.375
1.073 1.149 1.271 1.446 1.748 2.263 3.051 4.522
δ ) 0.3
δ ) 0.4
δ ) 0.5
η
T* ) 6.0
T* ) 2.5
T* ) 1.5
T* ) 8.0
T* ) 3.0
T* ) 2.0
T* ) 10.0
T* ) 4.0
T* ) 2.5
0.10 0.15 0.20 0.25 0.30 0.35 0.40 0.45
1.405 1.701 2.102 2.637 3.406 4.487 6.072 8.358
1.229 1.408 1.670 2.041 2.604 3.477 4.812 7.003
1.005 1.057 1.165 1.363 1.716 2.326 3.405 5.420
1.402 1.707 2.117 2.677 3.475 4.615 6.269 8.696
1.190 1.362 1.625 2.021 2.632 3.605 5.155 7.523
1.014 1.089 1.230 1.495 1.967 2.795 4.237 6.666
1.397 1.702 2.118 2.700 3.525 4.700 6.397 8.862
1.207 1.404 1.692 2.147 2.850 3.944 5.614 8.131
1.013 1.100 1.270 1.599 2.186 3.200 4.856 7.435
3724 J. Phys. Chem. B, Vol. 111, No. 14, 2007
Eu and Qin we examine the critical parameters (p/c , T/c , ηc) predicted by them. For this discussion we will not use the van Laar generalization8,9 for the definition of critical point, but the conventional one in which the first two density derivatives of pressure are set equal to zero. For example, if the virial equation of state, eq 58, is simply assumed applicable to the question, then we find there does not exist a physical critical point. However, if the equation of state is slightly altered by removing the denominator in the attractive term to the form
(1 + 2η + 3η2) p* ) - 12fδη(1 + η) ηT* (1 - η)2 f ) exp(1/T*) - 1
Figure 4. Comparison of theoretical results with MC simulations. δ ) 0.5. The meanings of symbols, the numerals, and curves are similar to Figure 1: 1 and open circles (O) for T* ) 10.0; 2 and triangle (4) for T* ) 4.0; and 3 and diamond (]) for T* ) 2.5.
regime. In the cited ref 12 it was shown that for an appropriate account of density dependence of critical isotherms in the neighborhood of the critical point there should be a fractional power term in the density expansion and that this fractional power term is directly related to the behavior of critical isotherms near the critical point: |p*/p/c - 1| ∼ |η/ηc - 1|4+R (R ) 0.3-0.5). The empirical GvdW parameters in the subcritical regime that can account for such a behavior are 3
A*(k) )
i 3 1+R (k) a(k) ∑ i (η - ηk) + Rna (η - ηk) |η - ηk| i)0
(R < 1, k ) l, V) (69) 3
B*(k) )
i 3 1+R (k) b(k) ∑ i (η - ηk) + βna (η - ηk) |η - ηk| i)0
(R < 1, k ) l, V) (70) (k) (k) where ηk (k ) l, V) are spinodal densities and a(k) i , bi , Rna , and (k) βna are parameters in the liquid (l) or vapor (V) branch, which are generally dependent on temperature and directly determined from the critical data, liquid-vapor coexistence curve, or spinodal curve, as discussed in refs 12 and 53. The perturbation theory results for the equation of state obtained, however, are regular with respect to η and T*. They, of course, can be expanded in series of (η - ηk), but the expansion does not give rise to such a fractional power (irrational) term, and neither is there a possibility for such a term to arise from the equations of state presented. As a matter of fact, it is not clear at present how such a fractional power term in the density expansion can arise in the integral equation approach if the PY closure is assumed. Therefore, it is evident that the PY closure should be improved in such a manner that the solution of the integral equation admits a singular behavior with regard to η and T*. This is a challenging question, and we hope to investigate it in the near future. Despite this weakness of the result, by gleaning the perturbation theory formula for the equation of state, or more simply, regarding the parameters * and δ as empirical parameters, it is possible to construct an empirical model for the equation of state to study thermodynamic properties and some aspects of fluid properties in the supercritical regime. To assess the qualitative aspect of capability of the equations of state obtained,
then there exists a critical point in the physically acceptable range:
ηc ) 0.163 T/c )
1
(
ln 1 +
0.715 δ
)
(fcδ ) 0.715)
p/c ) 0.062T/c Thus, for example, if we choose δ ) 0.2, then T/c ) 0.658, which suggests that T* ) 0.6 used for numerical comparison with MC simulation data in the numerical comparison section is indeed close to and below such a point in temperature, even if only the leading terms are retained in the attractive part. In this regard, it should be noted that the denominator of the attractive term does not alter the density dependence of the equation of state significantly at least in the range of relatively low η values in which ηc is located. On the other hand, if the compressibility equation of state is used, we find the critical point as follows:
ηc ) 0.204 fcδ ) 0.580 p/c ) 0.085T/c which implies T/c ) 0.477 for δ ) 0.1, 0.658 for δ ) 0.2, and 1.126 for δ ) 0.5. These values are off from the experimental values, but not completely off the range. Therefore, the analytic forms obtained for the equation of state can be made empirical provided the parameters are adjusted and, perhaps, another parameter or two are suitably added, so that the fluid properties near the critical point may be reproducible thereby. In any case, as mentioned earlier, these equations of state are not expected to properly account for the critical phenomena because of the absence of singular nonanalytic terms that causes them to be discontinuous in the subcritical regime and give correct critical exponents, and it is unrealistic to expect such a behavior from the PY integral equation, which is known to have a problem yielding converging solutions in the critical regime. To describe the subcritical behavior properly, it is necessary to assume the GvdW parameters, for example, (69) and (70). The rational power terms in these expressions for the singular GvdW parameters can be deduced from the perturbation theory results obtained in this work, but the irrational terms responsible for the discontinuity in the equation of state across the liquidvapor coexistence region cannot be obtained from the perturbation theory results in any way. They require a different treatment.
Square Well Potential Model
J. Phys. Chem. B, Vol. 111, No. 14, 2007 3725
We defer a study of how such an equation of state might arise in the integral equation approach to future work.
r‡3[1 + 2r‡3η* + 3(r‡3η*)2] p* ) η*T* (1 - r‡3η*)2
VI. Discussion and Concluding Remarks In this paper, we have solved the PY integral equation for a square well potential in the limit of small well width by treating it as a perturbation parameter and have obtained analytical forms of solutions. We have presented only the zeroth-order solutions for the cavity functions in closed forms and therewith calculated the GvdW parameters and also the virial and compressibility equations of state. The equations of states obtained are in analytic forms as functions of density and temperature. The repulsive part of the zeroth order solution for the cavity function is the well-known cavity function for the hard sphere fluid, and the corresponding repulsive part of the equation of state is the PY equations of state for hard sphere fluids via the virial and compressibility route. The lowest order cavity functions yield the attractive contribution to the first order in perturbation parameter. It, in fact, appears multiplied by fδ as if fδ is the perturbation parameter. The perturbation theory results are sufficiently accurate in the range of temperature T* J 0.6 for δ ) 0.05 (fδ ) 0.215), T* J 0.8 for δ ) 0.1 (fδ ) 0.249), T* J 1.5 for δ ) 0.3 (fδ ) 0.284), and T* J 2.5 for δ ) 0.5 (fδ ) 0.246) according to the test performed. These different cases have roughly similar values of fδ as shown, and it may suggest there is a value of fδ around which the PY approximation is no longer reliable. This does not imply that in the still lower temperature regime sufficiently away from the critical temperature the perturbation theory result may turn out equally unreliable; on the contrary, they might be as reliable as for temperatures higher than 0.6 or 0.8 or 1.2 or 2.5 tested. Further numerical tests should be made in the future in this regard. In a previous paper42 devoted to the behavior of the GvdW parameter B*, which is closely related to the excluded volume of the fluid, we have obtained an empirical formula for the temperature dependence of the diameter r‡σ of the excluded volume for the Lennard-Jones fluid, which in reduced units has the form
r‡ )
d0 (1 + d1xT*)1/6
(71)
where d1 is weakly dependent on density and the numerator d0 was found slightly larger than unity (see eq 11). We, however, find it possible to take d0 ) 1 and, to an approximation, d1 ) 0.76531. If the position variable r is scaled by r‡ instead of the contact diameter σ in the OZ-PY equations used under the assumption that the square well potential model has a temperature-dependent hard core, then the packing fraction for the excluded volume is given by
η* π η ) (r‡σ)3F ) 6 (1 + d1xT*)1/2
(72)
where η* ) (π/6)σ3F, the packing fraction of the hard sphere fluid of temperature-independent diameter σ. If this scaling is used and the perturbation theory is implemented, in the equation of state formulas the packing fraction is replaced by η defined by eq 72, and it would make the hard sphere part of the equation of state no longer athermal. For example, the virial equation of state for hard sphere fluids will have the form
(73)
Thus, if an empirical approach is taken for equation of state treating d1 as an adjustable parameter and perhaps using another exponent ϑ < 1/2 for temperature, that is, T*ϑ instead of xT* in eq 72, and if it is possible to add a suitable form for the contribution from the attractive potential, this kind of scaling argument may be used to invent a soft sphere equation of state that is not athermal. In connection with a utility of r‡, it must be noted that its use in the modified free volume theory formula for the self-diffusion coefficient for the Lennard-Jones fluid yields an excellent self-diffusion coefficient42 in comparison with experiment and simulation results. In concluding this work, the point we would like to make is that the analytic forms of equation of state provided by the perturbation theory for square well potentials supply useful semiempirical forms for the equation of state for study of thermodynamic properties of fluids away from the critical point, especially in the supercritical regime, and the considerable insight into the method of solution we have gained here may be used for study of similar solutions for complex fluids in the future. Last, there remains the challenge of improving the integral equations to be used in the present line of perturbation theory for thermodynamic properties. Acknowledgment. The present work was supported in part by the grants from the Natural Sciences and Engineering Research Council of Canada. References and Notes (1) van der Waals, J. D. Ph.D. Thesis, University of Leiden, 1873, A. W. Sijthoff, Leiden. (2) Berthelot, D. TraV. Mem. In. Poids Mes. 1907, 13, 113. (3) Beattie, J. A.; Bridgeman, O. C. J. Am. Chem. Soc. 1927, 49, 1665. (4) Redlich, O.; Kwong, J. N. S. Chem. ReV. (Washington, D.C.) 1949, 44, 233. (5) Peng, D. Y.; Robinson, D. B. Ind. Eng. Chem. Fundam. 1976, 15, 59. (6) Longuet-Higgins, H. C.; Widom, B. Mol. Phys. 1964, 8, 549. (7) See, for example: Guggenheim, E. A. Mol. Phys. 1965, 9, 199. (8) van Laar, J. J. Proc. Sci. Sec. Kon. ned. Akad. Wisensch. 1912, 14 II, 1091. (9) Baehr, H. D. Forsch. Geb. Ingenieurwes. 1963, 29, 143. (10) Eu, B. C.; Rah, K. Phys. ReV. E 2001, 63, 031303. (11) Eu, B. C. J. Chem. Phys. 2001, 114, 10899. (12) Rah, K.; Eu, B. C. J. Phys. Chem. B 2003, 107, 4382. (13) Eu, B. C. J. Phys. Chem. A 2006, 110, 831. (14) Eu, B. C. Transport Coefficients of Fluids; Springer: Heidelberg, 2006; Chapter 6. (15) Rah, K.; Eu, B. C. Phys. ReV. Lett. 1999, 83, 4566. (16) Rah, K.; Eu, B. C. Phys. ReV. E 1999, 60, 4105. (17) Rah, K.; Eu, B. C. J. Chem. Phys. 2001, 115, 9370. (18) Rah, K.; Eu, B. C. Phys. ReV. Lett. 2002, 88, 065901. (19) Volterra, V. Theory of Functionals and of Integro-Differential Equations; Dover: New York, 1959. (20) Alder, B. J.; Wainwright, T. E. J. Chem. Phys. 1957, 27, 1207. (21) Reiss, H. R.; Frisch, H. L. F.; Lebowitz, J. L. L. J. Chem. Phys. 1959, 31, 369. Helfand, E. H.; Frisch, H. L. F. J. Chem. Phys. 1961, 34, 1037. (22) Thiele, E. T. J. Chem. Phys. 1963, 39, 474. (23) Wertheim, M. Phys. ReV. Lett. 1963, 10, 321; J. Math. Phys. 1964, 5, 643. (24) Baxter, R. Aust. J. Phys. 1968, 21, 563. (25) Sastry, S.; Debenedetti, P. G.; Stillinger, F. H. Phys. ReV. E 1997, 56, 5533. (26) Debenedetti, P. G.; Stillinger, F. H.; Truskett, T. M.; Roberts, C. J. J. Phys. Chem. B 1999, 103, 7390. (27) Corti, D. S.; Debenedetti, P. G.; Sastry, S. Phys. ReV. E 1997, 55, 5522. (28) Barker, J. R.; Henderson, D. ReV. Mod. Phys. 1976, 48, 587. (29) Weeks Chandler, D.; Andersen, H. C. J. Chem. Phys. 1971, 54, 5237.
3726 J. Phys. Chem. B, Vol. 111, No. 14, 2007 (30) Titchmarsh, E. C. Theory of Fourier Integrals, 2nd ed.; Oxford: London, 1939. (31) Morse, P. M.; Feshbach, H. Methods of Theoretical Physics; McGraw-Hill: New York, 1953. (32) Noble, B. Methods Based on the Wiener-Hopf Technique; Pergamon: Oxford, 1958. (33) Duderstadt, J. J.; Martin, W. R. Transport Theory; Wiley: New York, 1979; Appendix A. (34) Hansen, J. P.; McDonald, I. R. Theory of Simple Liquids, 2nd ed.; Academic: San Diego, 1986. (35) Chapman, W. G.; Gubbins, K. E.; Joslin, C. G.; Gray, C. G. Fluid Phase Equilib. 1986, 29, 337; Fluid Phase Equilib. 1989, 52, 31. Mu¨ller, E. A.; Vega, L. F.; Gubbins, K. E. Int. J. Thermophys. 1995, 16, 705. (36) Dickman, R.; Hall, C. K. J. Chem. Phys. 1986, 85, 4108. Honnell, K. G.; Hall, C. K. J. Chem. Phys. 1989, 90, 1841. Yethiraj, A.; Hall, C. K. Mol. Phys. 1993 80, 469. (37) Wertheim, M. S. J. Stat. Phys. 1984, 35, 19, 35; J. Stat. Phys. 1986, 42, 459. (38) Flory, P. J. J. Chem. Phys. 1941, 9, 660. (39) Huggins, M. L. J. Chem. Phys. 1941, 9, 440. (40) Zwanzig, R. J. Chem. Phys. 1954, 22, 1420. (41) McQuarrie, D. Statistical Mechanics; Harper & Row: New York, 1976.
Eu and Qin (42) Laghaei, R.; Nasrabad, A. E.; Eu, B. C. J. Chem. Phys. 2006, 124, 154502/1-9. (43) Kirkwood, J. G.; Oppenheim, I. Chemical Thermodynamics; McGraw-Hill: New York, 1960. (44) Baxter, R. J. Chem. Phys. 1968, 49, 2770. (45) For papers discussing levitation effects, see, for example: Ghorai, P. Kr.; Yahonath, S.; Lynden-Bell, R. M. J. Phys. Chem. B 2005, 109, 8120. (46) Boltzmann, L. Lectures on Gas Theory; University of California: Berkeley, 1964 (translated by S. G. Brush). (47) Percus, J.; Yevick, G. J. Phys. ReV. 1958, 110, 1. (48) Rushbrooke, G. S.; Scoins, H. I. Proc. R. Soc. London 1953, A216, 203. (49) For a thermodynamic perturbation theory for a square well potential in this line of approach, see: Tang, Y.; Lu, B. C. Y. J. Chem. Phys. 1993, 99, 9828; J. Chem. Phys. 1994, 100, 3079; J. Chem. Phys. 1994, 100, 6665. (50) Carnahan, N. F.; Starling, K. E. J. Chem. Phys. 1969, 51, 635. (51) Allen, M. P.; Tildesley, D. J. Computer Simulations of Liquids; Clarendon: Oxford, 1987. (52) Binder, K.; Heermann, D. W. Monte Carlo Simulation in Statistical Physics; Springer: Heidelberg, 1988. (53) Eu, B. C. Unpublished results.
