Generalized mean spherical approximations of the hard-sphere system

Generalized mean spherical approximations of the hard-sphere system

Volume l14A, numbe~ 6 PHYSICS LETTERS 3 Maxch 1986 G E N E R A L I Z E D MEAN S P H E R I C A L A P P R O X I M A T I O N S OF T H E H A R D - S P ...

291KB Sizes 0 Downloads 59 Views

Volume l14A, numbe~ 6

PHYSICS LETTERS

3 Maxch 1986

G E N E R A L I Z E D MEAN S P H E R I C A L A P P R O X I M A T I O N S OF T H E H A R D - S P H E R E S Y S T E M Salvino C I C C A R I E L L O Dipartimento di Fisica "G. Galilei'" and Istituto Nazionale di Fisica Nucleare, sez. di Padova, Padua, Italy

Received 7 November 1985; accepted in revised form 19 December 1985

The known part of the direct correlation function outside the hard-sphere core is approximated by a sum of two yukawian contributions and the equations associated to the corresponding generalized mean spherical approximations are numerically solved. One finds that no physically acceptable solution exists when the dimensionless particle density exceeds the value 1.093. One concludes that, in this approximation, the fluid phase of the hard-sphere system cannot exist beyond the former density value.

L The hard-sphere (HS) model [1] is still the object of many theoretical and numerical investigations [2-6] due to its great importance for understanding the behaviour of real fluids. Among its features, undoubtedly the most fascinating one is the possibility of undergoing a liquid-solid phase transition in the density range 0.94-1,04 * 1. This phenomenon has been strongly suggested by the numerical calculations carried through by Adler and Wainwright [7] and by Wood and Jacobsen [8] in the middle fifties. Since then, many attempts have been made to obtain the former result in a more analytical way. On simple theoretical grounds in fact, this transition should be revealed by a mathematical singularity in the virial expansion of the pressure [9] at a particular density value PF. Unfortunately, only the first seven coefficients of the virial series are known exactly. In order to overcome this difficulty, people have assumed that the virial expansion is recovered by constructing Pad6 [2] or other approximants [3] with the known virial coefficients. The strength of this assumption is evident. Other people [10,6] have looked for a possible mechanical instability of the system. More explicitly one determines the trajectory of the complex root z(p) of the equation ,2 , I We are using the dimensionlessparticle density (---Nda/Vvalues. Moreover,from now on, the unit of length will be

p )

taken equal to the particle diameter d.

0.375-9601/86/$ 03.50 © Elsevier Science Publishers B.V. (North-Holland Physics Publishing Division)

1 - ~ ( z , p) = O.

(1)

If the z(p) trajectory meets the real axis at a density value, say PMI, smaller than the close-packing value, the system is said to be mechanically unstable and a phase transition must take place at a density not larger than PMi *a . Besides these two ways, we have another method for attempting to locate the density where the existence of the fluid phase breaks down. In the case of HS in fact, the virial route to thermodynamics relates the pressure to the contact value of the radial distribution function (rdf): g(1 +, p) and this quantity, in turn, is equal to the discontinuity of the direct correlation function (dcf) at r = 1. According to Mayer's analysis, the onset of freezing should be revealed by the appearance of a singularity in (c(1 +, p) c ( 1 - , p)), too. Of course the search of a singularity is meaningful only for the dcf approximations which are subseries of the exact cluster expansion. These dcf expressions can be immediately obtained from any ,2 In oq. (1) ~(z, p) denotes the analytic continuation over the complex z-plane of the Fourier transform ~(k, p) (FT) of the direct correlation function. In this paper we shall adopt the notation of ref. [6] which will be quoted as I in the following. .a Actually, present experimental results suggest that this density must be smaller than PMI because the peaks of the structure function do not show any dramatic increase as the liquid approaches to the freezing point [ 11 ]. 313

Volume 114A, number 6

PHYSICS LETTERS

integral equation which has a clear cluster meaning. It is important however that the resulting dcf's have an expression as close as possible to an algebraic one. This happens for the class of the so-called generalized mean spherical approximations (GMSA) [12,13]. In that case the r-dependence of the dcf is known [see eq. (2)], while the density dependence is present only in some unknown numerical coefficients, i.e.: a, b, u 1, .... oN in eq. (2). Moreover the latter ones must obey a system of algebraic equations. Therefore if one finds that, by letting the density increase, the equations do not admit any more real solutions beyond a particular density value (OB), one concludes that the considered dcf approximation is singular at PB and that the pure fluid phase cannot exist beyond OB. Of course, in order for the result to be considered reliable, the resulting approximation must yield a thermodynamical and structural picture of the system which looks globally satisfactory. 2. We can now sketch how one can practically carry through the aforesaid program. Examples of integral equations with an evident cluster meaning are the ones discussed in section 2 of I. They correspond to applying the iteration procedure given by eqs. 0-2.4, 2.5) and by assuming that, except for a discrete set of values of n, all the C ~ ( r ) are null. [ C ~ ( r ) denotes the coefficient of pn in the density expansion of Cext(r, p).] These equations will be referred to as generalized Percus-Yevick equations*4. We recall now that Cext(r, p) is exactly known only up to terms o(p3). Besides, C(e2~(r) and C(e3x)t(r)have their supports in (1, x/~) [15] and (1, 1.8) [16], respectively, while their shapes are such that they can be well approximated by a sum of yukawian terms. With this approximation, the generalized PY equations, obtained from the aforesaid two known expressions, become equivalent to the relevant GMSAs. These approximations guarantee that Cint(r, p) has the following functional dependence [13] -Cint(r, p) = a(1 + 0.5~r 3) + br 2 N

+ E viV(zi, gi, oi, r), i=I

3 March 1986

N Cext(r, p) = ~ Ki(P) exp[-zi(r - 1)]/r. i=1

The definition of V( ) can be immediately obtained by recalling that eq. (2) corresponds to eq. (la) of ref. [13]. a, b and o I ...uN are unknown quantities which must be determined by solving the system of equations (2) of ref. [13]. The analytical expressions of eqs. (2) and (3) yield simple analytical expressions for the FTs of the dcf's. The resulting approximations of d(k, p) can be analytically continued in the complex z-plane and one can look for the first complex zero of eq. (1) along the same way of I. Moreover the use of the OrnsteinZernike equation yields the FT h(k, p) of the total correlation h(r, p). Finally, by inverting h(k, p), one can determine h(r, p) and obtain a further test on the resulting accuracy of the considered GMSAs. 3. After having expounded the general line of the work, we can now briefly illustrate the results, leaving to a future paper a more detailed report. We have approximated C(e2)(r) both by one and by two yukawian terms, while C~)t(r ) has been approximated by only one term. The values of parameters (if} and {X} appearing in the approximations

Ni C(e/x)t~ Ei -j~-I ki/j exp[-Xi/J (r - 1)]/r,

i = 2, 3 ,(4)

have been determined by requiring that E i (1 +) = C(e/) (1+) and that the L 2 norm o f E i - C(iex)tbe as small as possible. The parameter values are reported in table 1. We have considered the GMSAs relevant to the following cases: TF1, FT1, TF2 and TT2. In order to

Table 1 Values of (/( } and {k} which yield the best norm value (second column) for the difference function relevant to the approximations shown on the left.

(2)

L2

g

k 5.998 9.544 10.645 8.607

once one has assumed that

TFI

0.518 X 10 -3

,4 In fact their solutions fulfill exactly the core condition, are exact up to o(p N) when o n e u s e s the exact C(enx) t with n
TF2

0.621 X 10 .-4

0.1621 1.0307 -0.8686

FT1

0.140 X 10 -1

0.382

314

(3)

Volume l14A, number 6

PHYSICS LETTERS

understand the meaning of these symbols we recall ^2,'~(2)(r~ that the known part of Cext(r, p) is ~, ~extv J + p3r,(3). (rI Thus in TF2, for instance, the first letter ~ ext,k /. (T) means that the first contribution (C(e2x))has been accounted for, the second letter (F) means that C(3ex)t has been neglected, while 2 indicates that Cext(r, p) has globally been approximated with two Yukawa terms. We have numerically solved eqs. (2) of ref. [14] relevant to all the four aforesaid cases. We have determined the roots a, b, Ol, ..., oN in the density range [0;1.5], but for the TT2 approximation, where we have not been able to go beyond PB = 1.092. One finds that the thermodynamical consistency improved as we go along the approximation chain TF 1, TF2, FTI, TT2, as is clearly shown by the contact values reported in table 2 *s. We can now report the results relevant to case TT2. Fig. 1 shows the trajectories of the first complex zero of the eqs. (1) associated to the analyzed GMSAs. , s The corresponding radial distribution functions, obtained via a fast FT algorithm, once they are compared with the MC ones of ref. [17], show a similar improvement only in the contact region. By contrast, the agreement become worse in the regions where maxima and minima occur. The amplitudes of these in fact, are larger than the corresponding MSA ones and increase along the aforesaid chain. We remark that all the results relevant to case TF1 are close to the TF2 one as one would expect since C(e~t(r)turns out to be already well approximated by one yukawian term. We have calculated all these quantities in order to be humanely confident on the absence of computer progzam mistakes.

3 March 1986

9= 1.05

)=1,10

0.30 Im(z) 0.20F

\ \~Q=tzo

0.10

72

7.4

7.6 7.8 Re (z)

8.0

8.2

8.4

Fig. 1. The trajectories of the first zero of eq. (1) for approximations: MSA (continuous line), TF2 (broken line), FTI (broken bold line) and TT2 (broken dotted line) are reported. The trajectory relevant to the case TF1 cannot be distinguished from the TF2 one.

Table 2 The contact values relevant to the GMSAs considered in the paper are reported in columns 4-7, while for comparison columns 2 and 3 show the MC and the MSA [17] values. Note the small difference between the TF1 and the TF2 approximations. MC

MSA

TF1

TF2

FT1

TT2

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9

-

1.143

1.144

1.144

1.143

1.145

1.339 1.548 1.811 2.116 2.561 3.157 3.971 5.147

1.0

-

1.313 1.518 1.768 2.075 2.460 2.949 3.581 4.419 5.559

1.319 1.532 1.794 2.124 2.546 3.093 3.805 4.716 5.881

1.319 1.531 1.792 2.117 2.531 3.068 3.777 4.722 5.964

1.316 1.527 1.789 2.120 2.545 3.103 3.856 4.888 6.286

1.322 1.541 1.817 2.171 2.634 3.247 4.036 5.051 6.686

315

Volume 114A, number 6

20

l

I

PHYSICS LETTERS I

I

I

I

I

10

3 March 1986

On the other hand if one recalls that a, b, o 1, and v2 can be expanded in a power series in p, the former result implies that at PB the expansions have an analytical singularity which from eq. (2) and by the argument expounded in section 1 will be present also in the virial expansion of the pressure* 7. Therefore one can conclude that, at least in the TT2 approximation, the fluid phase of the HS system cannot exist beyond PB. Moreover, both the fact that at PB the HS system is still rather far from being mechanically unstable and the other numerical results indicate that approximation TT2 looks satisfactory enough for making the former result reliable. ,7 I thank a referee for a remark related to this point.

-10

References

./ -20 0.0

*

I

OA

I

I

0.8

I

I

1

1.2

0

Fig. 2. Physical {'~1.....Y~4}and unphysical {'~1..... "~t} roots of the TT2 GMSA equation. It appears evident that they will meet at PB = 1.092 -+0.001. The Xs and the Xs are related to a~ b, oI and 02 in the following way: X1 = -(aMS A - a) + 9, X 2 = -(bMs A - b)/6~ + 18,-~a = 5 X 104Vl/(24~p2k211i(2/l e_xP0h.2/,) - I and X4=5 X lOgv2](24~p31(3/lh3/1 exp h3/1)

One sees that in contrast to I, where the core condition was not fulfilled, the mechanical stability is recovered for approximations TF1, TF2 and FT1. By contrast, the trajectory relevant to the TT2 case ends at PB = 1.092. The reason of this fact is clear from fig. 2 where we have drawn the trajectory of the roots of the corresponding GMSA equations. One sees that at PB the physically acceptable root meets with one of the unphysical ones, so that beyond the former density value the root will become complex .6 and thus physically unacceptable. This mechanism is essentially the one noticed by Baxter [20] for the first time in the adhesive HS model. ,6 We recall that the system of GMSA equation can be written in a form which makes its algebraic nature evident [18,19]. 316

[ 1] J.A. Barker and D. Henderson, Rev. Mod. Phys. 48 (1976) 587. [2] G.A. Barker, G. Guttierrez and M. deLlano, Ann. Phys. (NY) 153 (1984) 283. [3] A. Baram and M. Luban, J. Phys. C12 (1979) L659. [4] B.B. Deo and A.C. Naik, Phys. Rev. A28 (1983) 1700. [5] G.A. Martynov and G.N. Sarkisov, Mol. Phys. 49 (1983) 1495. [6] S. Ciecariello and D. Gazzillo, Mol. Phys. 54 (1985) 863. [7] B.J. Aider and E. Wainwright, J. Chem. Phys. 27 (1957) 208. [8] W. Wood and J.D. Jacobson, J. Chem. Phys. 27 (1957) 1207. [9] J.E. Mayer, in: Encyclopedia of physics, Vol. XlI, ed. S. Fltigge (Springer, Berlin, 1958) pp. 73-204. [10] K. Hixoike, J. Phys. Soc. Japan 13 (1958) 1497. [11] A. Sj61anderand L.A. Turski, J. Phys. Cll (1978) 1973. [12] E. Waisman, Mol. Phys. 25 (1973) 45. [13] J.S. Hoye, G. Stell and E. Waisman, Mol. Phys. 32 (1976) 209. [14] M.S. Wertheim, Phys. Rev. Lett. 10 (1963) 321; E. Thiele, J. Chem. Phys. 38 (1963) 1959. [15] B. Nijboer and L. van Hove, Phys. Rev. 85 (1952) 777. [16] F.H. Ree, R.N. Keeler and S.L. McCarthy, J. Chem. Phys. 44 (1966) 3407. [17] J.A. Barker and D. Henderson, Mol. Phys. 45 (1971) 397. [18] J. Konior and C. Jedrzejek, Mol. Phys. 48 (1983) 219. [19] P.T. Cummings and E.R. Smith, Mol. Phys. 38 (1979) 997. [20] P.J. Baxter, J. Chem. Phys. 49 (1971) 2770.