Cohen, E. G . D . Rethmeier, B.C. !958
Physica XXIV 959-969
CELL-CLUSTER
THEORY
OF THE LIQUID
A TWO DIMENSIONAL
STATE. V
FLUID OF HARD SPHERES
b y E. G. D. COHEN and B. C. R E T H M E I E R (Instituut voor theoretisehe Physiea, Universiteit van Amsterdam, Nederland)
Synopsis The cell-cluster t h e o r y of the liquid state, developed in two previous papers i) 2), is applied here to a two dimensional fluid of hard spheres. The corrections introduced by considering cell-clusters of two cells are investigated. The double cell p a r t i t i o n function is calculated in an a p p r o x i m a t e - and an e x a c t w a y and the influences on e q u a t i o n of state and e n t r o p y are m u t u a l l y compared. Finally the two results of the cell-cluster t h e o r y are c o m p a r e d w i t h the Monte Carlo calculations for a two dimensional fluid of h a r d spheres by M e t r o p o l i s , R o s e n b l u t h , Rosenbluth, T e l l e r and T e l l e r 6).
§ I. I n t r o d u c t i o n . In two previous papers (referred to as I 1) and 112)) a systematic method has been given to extend the Lennard-Jones-Devonshire single cell theory of the liquid state 3). Thereto instead of only one cell, clusters of in general l cells were considered in which l molecules are supposed to move in a potential field composed of 1) their mutual interaction and 2) the interaction with the other N -- l molecules placed at the centres of their cells. Indicating the partition function of the molecules in a cluster of l cells ~...~ as Q'~...A, v-functions can be defined in terms of these partition functions b y the relations:
Q'~...~ = Ed~v v~...~ v,...~ ... %..~
(1)
where the summation extends over all possible divisions of the cluster ~... 2 into clusters of equal and smaller size. From these relations the v-functions can be solved to give: v l = v~ = Q~' = QI' v2 = v,~fl = Q~fl' .
.
.
,
°
.
.
.
.
.
.
.
.
.
.
. * *
. . . .
.°
.
(2a)
Q~'Qfl' = Q~' .
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
Q1 '2 .
.
.
.
(2b)
. . .
as has been explained fully in I and II. If we now restrict ourselves to clusters of one and two cells only, we can write the total partition function ZN in terms of the functions vl ata.d va in the following form: Z N = Z N (L'J', [(1 - - #,)z/2/(1 - - Z#,)] N --
9 5 9
--
(3)
960
E.G.
D, C O H E N A N D B. C. R E T H M E I E R
with =
2zw2 + 1 - - a / 4 ( z - - 1)w2 + 1 2(z2w2 + 1)
(3a)
where ~.~zN gives the n u m b e r of double cells t h a t contribute to the partition function and where w2 = vs/vl 2. ZN ~L'J'~ is the single cell partition function as used b y L e n n a r d - J o n e s and D e v o n s h i r e . Using for the Helmholtz' free energy F = -- kT In Z2v and, using the fact t h a t for a s y s t e m of hard spheres wg. does not depend on T, we find for the e n t r o p y and the pressure: S = S (L.J.)
2V N k In {(1 -- $)~9'/(1 - - za~)}
z(1 + (z- 2),~} c1~ ( awg.'~ p ___ #(L.J., -F kT 2(1 -- a)(1 -- za) dw2 \--~--v/~'
(4) (5)
respectively. In a foregoing paper 4) we investigated the influence of the correction term due to double cells for a hypothetical liquid composed of hard elastic spheres. In t h a t case the intermolecular potential is the unit step function: exp - - fig(r12) = ~'(r12 - - a) where 8 ( r 1 9 . - a) = 0 for r12 < a #
~ ' ( r l 2 - - a ) = 1 for r 1 2 ~ a so t h a t the partition function for single- and double cells can be written as:
QI' = ./singlecell drl
Q2' =-- ½ffaouble cell •(r12 -- a) d r l dr2
(6) (7)
respectively. Here the integrations are carried out over the volume left for the molecule(s) in the cluster between the collision spheres of the surrounding molecules. For the calculations of the single cell- or free volume contribution ZN ~L'd'~, the sphericalized single cell was used, which gives: QI' =
-
(8)
where aS = y V / N = yv is the relation between the distance of neighbouring cell centers a and the volume per molecule v = V/N. ~ depends on the structure of the lattice of the cell centres. For the calculation of the correction due to double cells,the complicated integration region of the double cell integral Qs' was approximated by the
A TWO D I M E N S I O N A L F L U I D OF HARD S P H E R E S
961
inscribed parabmepipedand a similar approximation was used for the single cell integral QI' in this case. This approximation of the real cell volume b y a parallelepiped is especially bad at high densities, as is illustrated in figure 1b. Still one can hope that w~ -- (Q2'/Q1 '2) - I is given essentially in a correct way. Calculations made for a simple cubic lattice (? = I) and a cubic face centered lattice ( ? - I/~/2) were compared with the results obtained with the Monte Carlo calculations b y R o s e n b l u t h and R o s e n b l u t h 5) for a liquid of hard spheres. The results showed that for the equation of state the corrections of the free volume theory calculated on the basis of the cellcluster theory were qualitatively in agreement with the Monte Carlo results, but for a shift of the whole of the cell-cluster curve in the direction of the positive v*-axis with respect to the Monte Carlo results. As was pointed out before 4), this could be due to: (a) the approximations introduced in the evaluation of the single- and double cell partition functions; (b) our restriction to double cells; (c) the possibility of vacant sites (holes) in the lattice of the cell-centres. The effect of these vacant sites Will be considered in a later paper. In this paper we investigate the effect of the approximations listed under (a), b y comparing an approximate with an exact evaluation of the single- and double cell partition functions of a two dimensional fluid of hard spheres. In the next section we give the results for an approximate treatment of the two dimensional fluid, which is similar to the one given in the previous paper for the three dimensional fluid. In § 3 a method is described to evaluate the single and double cell integrals exactly and § 4 compares the results obtained with this method with those of § 2.
§ 2. Approximate evaluation o / t h e cluster integrals/or a two dimensional hard sphere fluid. The previous calculations for a three dimensional fluid have shown that only a cubic face centered lattice as a basic structure for the cell-centres leads to agreement with the Monte Carlo results. Correspondingly we chose for the two dimensional lattice a trigonal lattice as a basis. The geometry of the approximate single and double cells are indicated in figure l a and in l b for four characteristic values of the reduced volume. In the rectangular approximation the double cell integral (7) can be written as : Q2' = (ll2!) foa' foa' f ~ foa~ ~ (r12 -- a) dx2 dy2 dXl dyt where, introducing y = (a/a)--
1 =
a*
--
(9)
l:
dl* = dlla = (a/a) ~/5 -- 2 = 1.732 (y + l) -- 2 d2* = d~./a = (3a/a) -- 2 = 3y + 1
(lOa) (lOb)
962
E. G. D. COHEN AND B. C. R E T H M E I E R
T h e algebraic i n t e g r a t i o n results are g i v e n in t h e a p p e n d i x . F o r t h e p u r p o s e of e v a l u a t i n g {)2' - - {)1'2 it seems to be m o s t consistent to e v a l u a t e Q I ' b y a p p r o x i m a t i n g t h e single cell also b y a r e c t a n g l e w i t h sides 2a - - 2g a n d dl = a ~ / 3 - - 2~. So {)1'* = {)1'/~ 2 = (2a* -- 2 ) ( a * ~ / 3 - - 2) = 2y {1.732(y + 1) - - 2}
(11)
I n t a b l e I one finds t h e v a l u e s of w2 ca) = v2/vl 2 a n d dw2~")/dy for t h e a p p r o x i m a t e cell as a f u n c t i o n of y a n d t h e c o r r e s p o n d i n g values of t h e r e d u c e d t o t a l a r e a A * = A / N g 8 = (~/3/2)(y + 1)2. T h e s u p e r s c r i p t (a) i n d i c a t e s t h e a p p r o x i m a t e n a t u r e of this e v a l u a t i o n of w2. No values of w2 ca) c a n b e g i v e n for v a l u e s of y smaller t h a n y = (2/~/3) - - 1 = 0.1547, b e c a u s e for this v a l u e of y t h e double cell d e g e n e r a t e s into a line in t h e r e c t a n g u l a r a p p r o x i m a t i o n (fig. lb). T h e v a l u e of 1/8 f o u n d for w2 (~) in this limiting case, is t h e w2-value for a one-dimensional fluid. tt
i
T
',
Io
....L
/
,,
1 ..~
/6
/
I
I
o
o
II I
/
I
,L -r
.~ "
/o
I
I
l ~
I'I l o,!
~ -
o
I
I
£ T
I I
o "
\o
\
Fig. le. The figure shows part of a two dimensional trigonal lattice, with lattice constant a. The drawn circles are the collision circles with radius a(a is the diameter of the hard spheres). They enclose the double cell in which one finds the two lattice points e and ~. To get the two single cells we draw the collision spheres round e and ~ (broken cii'cles). The - - . - - . lines indicate the boundaries of the approximate single and double cell, the latter with length d = 3 e - 20 (single cell 2 a - 20) and height d = a ~/~ -- 2a (same for single cell). The hexagon (broken lines) represents the are.a available for each molecule in the limit of negligible diameter: a--~ 0.
§ 3. Exact evaluation o / t h e cluster integrals/or a two dimensional fluid o/ hard spheres. T h e g e o m e t r y of t h e t r u e single- a n d double cells are s h o w n in figure l a. T h e a r e a of t h e e x a c t single cell c a n be e v a l u a t e d to be:
gl'*:
3
3
(y + 1)2 ~/~ _ ~_ (y + 1) V'4 - - (y + 1)~ - - b arc sin y +2 1 + ~ (12)
9~3
A TWO DIMENSIONAL FLUID OF HARD SPHERES
~,,.'~JJFI I]1 I ] I II~...,,~ •L..,~--I--I~ - ~ . -I- ..I-.-I- ..1-4 -~- ~ " p
\I
II
I
~
L~,~:I I I I I I|I I I I
.II~:.-.,..,,,,,,.
~~.j,~l'ly,~,]I I (I11~ lIl lI l lI l Il l Il lIII, , i ,I JI
~~J
L
IIIII/I
Jv
I I I I II
1.
II
y=O.05 Fig. lb. I l l u s t r a t i o n of t h e e x a c t e v a l u a t i o n of w~ for f o u r v a l u e s of y, as e x p l a i n e d in t h e t e x t .
~)64
E.G.D.
C O H E N A N D B. C. R E T H M E I E R
for 0 _< y _< ~/3 - - 1 = 0.732. F o r y = 0.732 the collision spheres of the nearest neighbour molecules just intersect each other at the c o m e r points of the hexagon (see fig. 1a), which determines the single cell. For larger values of y the molecule is n o t really confined to m o v e in an area of A / N as it should be. Thus we confined ourselves to values of y in the range given above. The calculation of ws was performed b y a graphical m e t h o d b y plotting the double cell, as pictured in fig. 1b, for several values of y = a* -- 1. The fourfold integration over the positions of the molecules I and 2 was t h e n carried out in three steps: First we drew the collision sphere a r o u n d the centre (Xl, Yl) of the molecule 1 and e v a l u a t e d with a planimeter the shaded area as indicated in fig. l b, counting it positive or negative respectively as indicated in the figure. This shaded area was determined for a set of 100 values of Xl, yl arranged in a network, laid over one q u a r t e r of the double cell (fig. l b). T h e resulting f u n c t i o n / ( X l , Yl) was t h e n integrated for each value of Xl from Yl = 0 up to yl = yb(Xl) (y~ is the m a x i m u m value of y for a given xl) and finally this result was integrated from Xl = 0 up to the m a x i m u m value of Xl occurring in the double cell. I t is clear t h a t the result of the integration gives the value of ¼(Qg.' - A 2Ql') where A2 is the total area of the double cell. The values of wg. = (Q2' - QI's)/ Q1 '~ can t h e n be found i m m e d i a t e l y b y a careful d e t e r m i n a t i o n of A ~ and QI'. The values of w2 were determined in this w a y for values of y = 0.3 up to 0.7 with intervals of 0.1. W h e n y < ( 2 / ~ / 3 ) - - 1 = 0.1547, the double cell is cut into two parts because of the overlapping collision spheres, as indicated in fig. lb. In this region two other determinations of we have been m a d e for values of y = 0.1547 a n d y = 0.05, which, however, are less accurate t h a n the previous ones. The limiting values of Q2' and QI' for y -+ 0 can be determined analytically to be Q2'* = (217/18)y 4 and QI'* = 2~/3y 9", leading to a limiting value we = 1/216 = 0.00463 for y = 0. I t is interesting to note t h a t Qg.'* does not approach Q1 '.2 in the limit y = 0 because of the occurrence of small trian~,mlar regions (see fig. l b), which continue to contrib u t e to Q2'* even in the limit y = 0. I t is difficult to get a reliable estimate of the accuracy of the calculations. The accuracy of the integrals for y > 0.1547 is probably about 5~oo. For y < 0.1547 the accuracy of the calculations becomes m u c h poorer. The results are summarized in table II. Here one finds the values for wg.~*~and dw~t~/dy as a function of A* = (~/3/2)(y + 1) 8. The derivative dwz~tl/dy was determined graphically from the curve giving w2 ~t~ as a function of y. The supercript (t) indicates t h a t the values refer to a result obtained w i t h the true cell.
§ 4. Results and discussion. The results obtained for wg.~J and wg.~t~ given in the tables I and II, are plotted in figure 2 as a function of A*~ = (A/Na~') ~-. It is clear, t h a t the agreement between the two curves is not so good,
A TWO DIMENSIONAL FLUID OF HARD SPHERES
96"5
although the rectangular approximation gives for larger area's a value for w~.ca~ of the correct order of magnitude and of the correct general trend. The two curves differ in particular in the region of small area's. The area for o
~ ~ ( t )
O:
s S~ • ~'r.~ ~
o Fig. 2:w2
=
'
T
12
1£
'~
Q 2 ' / Q z '2 -
~4
"
li6
o:~
fromthe exact(t)
1
,
o'.e
~ll
'
1~'
--.y
and the approximateevaluation(a).
-
The arrow indicates the value for which the molecules can pass each other for the first time. Below y = 0.3, w2 m is dotted because of poorer accuracy.
~S
T ~3
(12
~ ~ { t )
it
iI iI
0.1
/ / iI ii
°~1'°'
.
'o:2 1~2
F i g . 3. z I S / N k = ( S - - S t L . J . * ) / N k
•
~
o.,~
1/
o:~
from the exact
'
1,,
(t) - a n d
o;e
.
1,8
--~y
to
approximate (a) evaluation.
which the two molecules can pass each other in the double cell for the first time is much lower for the real cell A~,t* = 1.84 (y = 0.46) than for the rectangular cell A~,a* = 2.60 (y = 0.73). The increase of w~~t~ for A* > 1.84
966
E. G, D. C O H E N
AND
B. C. R E T H M E I E R
is much steeper than that of w ~ (a) for the corresponding values A* > 2.60. A very similar behaviour is shov~n b y the entropy correction A S = = (S -- s(r"J'))/Nk as a function of A*½, which can be calculated using the values of w2 (al and wg.(t) as given in the tables I and n respectively and (4). The numerical values are given in the fifth column of the tables I and n a n d t h e y are plotted in figure 3. TABLE I
0.15 0.20 0.30 0.40 0.50 0.60 0.70 0.80 0.90 1.00
0.125 0.127 0.138 0.154 0.173 0.194 0,219 0.255 0.293 0.323
0 0.082 0.142 0.172 0.198 0.229 "0.287 0.402 0.345 0.251
1.155 1.248 1.464 1.698 1.949 2.217 2.503 2.806 3.126 3.463
0.243 0.246 0.262 0.282 0.302 0.327 0.354 0.388
0.422 0.446
0 0.067 0.120 0.149 0.172 0.199 0.246 0.333 0.276 0.198
7.464 .6.067 4.454 3.649 3.172 2.866 2.674 2.583 2.389 2.198
7.533 6.140 4.537 3.745 3.282 2.995 2.828 2.771
TABLE II
I ws(° 0 0.05 0.1547 0.30 0.40 0.50 0.60 0.70
[
du"(~)/dY
0.00463
(0.05)
(0.0089) (0.0311
(0.16) (0.30)
0.0711 0.0818 0. l 108 0.1763 0.2473
0.219 0.088 0.530 0.715 0.655
I A* 0.866 0.955 1.155 1.464 1.698 1.949 2.217 2.503
I ASU)/NI¢ [ APA(*)INkT [ Pt°AINkT I P~U)AINkT 0.0135
(0.026)
(0.080) 0.160 0.179 0.224 0.308 0.381
(0.07) (0.23) (0.39) 0.203 0.I12
oo
oo
(21.23)
(7.85)
0.569 0.686 0.517
4.536 3.612 3.569 3.353 2.946
(21.291 (7.92) 4,620 3.708
3.679 3.482
3.099
The corrections on the equation o~ state A p A / N k T = (p -- p ( L . J . ) ) A / N k T for the two cases calculated from equation (5) are also given in table I and n and plotted in figure 4. The shape of the two curves is very different. In the same figure the values for A p A / N k T are plotted derived from equation of state calculations with the Monte Carlo method by M e t r o p o l i s , R o s e n b l u t h , R o s e n b l u t h , T e l l e r and T e l l e r for a two-dimensional hard sphere fluid 6) by subtracting from their results the values of the single cell approximation: p ( L . J . ) A / N k T = (y + 1)/y (14)
calculated with a spherical single cell. The pressure determined from the exact expression (12) for the single cell is pt(L'J')A NkT --
0 In QI' (y + l) 0 In A -- 2Ql'* [3%/3(y + 1) - 3%/4 -
(y + 1)9'] (151
valid for y < 0.8. The values obtained with (14) and (15) show little difference.
A TWO DIMENSIONAL FLUID OF HARD SPHERES
967,
Fig. 4 shows, that whereas the pressure corrections on the single cell theory determined from the cell-cluster treatment appear to be positive in both calculations, the pressure obtained in the Monte Carlo calculations is 1.2
T 1.0 O.8 •
0.6
",
it)
0.40.2- /
d3
c
6
-o.~, -11.4 -I1E
Fig. 4. A p A / N k T = (p -- p~L.J.OA/NkT f r o m e x a c t (t) and a p p r o x i m a t e (a) e v a l u a t i o n c o m p a r e d w i t h a curve calculated f r o m the M o n t e Carlo values b y s u b t r a c t i n g f r o m p~M.O.JA/NkT the pcL.J.)A/NkT = (y + !)/y.
T 8
7 6
~..
5
3 2 1 0
• 1.o I ,
nl
1.1 I,
02
1~,
03
1~I,
Q~
1.~ ,I
QS
~,~ aS
~i6 (Z7
1.7 li8 Qe'
-.A *~' 'LO--.y
Fig. 5. p A / N k T curves as calculated by adding J p ' t ' A / N k T and .4p~a).4/NkT resp. to p~L.J.).4/NkT = (y + 1)/y.
smaller than that obtained in the single cell approximation. The authors have been unable to find an explanation for this. It should be mentioned of course that, although the evaluation of wg(A*) in § 4 is exact, the solution
968
E. G. D. C O H E N A N D B. C. R E T H M E I E R
of the combinatorial problem as given in the equations (3a) and (5) involves approximations. It m a y be, however, "that the origin of the discrepancy is to be found in the Monte Carlo results for which the estimated possible error in the absolute values of p A / N k T becomes of the order of 0.4 in the region of A* = 1. In figure 5 finally the values of p A / N k T are plotted. They are tabulated in the tables I and II; also are given there for comparison the values of p A / N k T obtained b y adding A p m A / N k T and A p c ~ M / N k T to both plL..1.M/NkT and ptlL.J.M/NkT. The authors are much indebted to Professor Dr. J. de B o e r for many stimulating discussions. Received 12-7-58. APPENDIX
Calculation o/ Q2" in the rectangular approximation *) We shall evaluate the fourfold integral Q2' = (1/2!) foal foa*foal foa* .~ (r12 - a) dx~ dy2 dXl dyl
(A I)
where r12 2 =
(x2 -
xl) 2 +
(y~ -
yl) ~
Introducing a set of relative coordinates: X
=
xg. - - X l
( A 2)
Y = y2 -- Yl • and writing x, y for xl, Yl respectively we can write (A 1) in the form :
Q2' = ½foa'fa-~ -~ K ( Y ) d Y dy
(A3a)
where
K(Y)
_- - J Ora, ras-x @(R -- a) d X dx d--x
(A 3b)
with R2 =
X ~. +
y2
In each of the above double integrals we change the order of integration. This leads to: Q2' = 2 foa*foa' N ( X , Y) o~ (R -- a) d Y d X (A 4) with N ( X , Y) = (d2 -- X)(dl -- Y) (A 5a) or
N(R, v~) = did2 -- (dl cos 0 + d2 sin v~) R + R 2 cos ~9sin v~
(A 5b)
Here the second expression is given in terms.of spherical polar coordinates. Handling the integral (A 4), two cases have to be distinguished: *) This calculation is due to Professor Z. W. S a l s b u r g .
969
A TWO DIMENSIONAL FLUID OF HARD SPHERES
C a s e I dl > a or y > 0 . 7 3 2 The integral can be c o m p u t e d from:
Q2'cz~ = ½dlZd22 -- 2 f ~ 12 dO f ~ N(R, ~9) R dR
(A 6)
(Q2 q]~ = ½ d12d22 when a = 0) as is illustrated in fig. 6. One obtains:
Qg.'*~I) = ½dl*2d2 .2
O"
I
d=
-
½~dl*d2* + ~(dl* + d2*) -- ¼
-
,X
¢~ d 2
'
(A 7)
~X
(a) Cb) Fig. 6. R e g i o n s of i n t e g r a t i o n in t h e t w o d i f f e r e n t cases for which Q2' h a s t o be e v a l u a t e d :
(a) case I: a < dl, (b) case II: dx< o < d~ C a s e II de > ~ > d l or 0 . 1 5 4 7 < y < The integral is represented b y :
0.732
Q9.,,II, = Q~,,I, + 2f~l~ d0f•xeoseo o N(R, ~9)R dR
(A 8)
as is shown in fig. 6. Here 01 = arc sin dl/a The result of the integration is:
Qg.'*IH~ = ½dl*~'d~.2 + }d2* -- ~ d1.4 -- dl*d2* arc sin dl* + --
}d2*(½d1.2 + 1) ~/1 -- d1.2 + {d1.2 (A 9)
The expressions for dl* = dl/a and de* = d2/a as functions of y are given b y (10a, b). Received 12-7-58. REFERENCES I) 2) 3) 4) 5) 5)
De Boer, J., Proc. conf. theor. Phys. Japan, 1953, 508; Physica 20 (1954) 655. Cohen, E. G. D., De Boer, J. and S a l s b u r g , Z. W., Physica 21 (1955) 137. L e n n a r d - J o n e s , J. E. and D e v o n s h i r e , A. F., Proc. roy. Soc. A 163 (1936) 63; 164 (1938) I. S a l s b u r g , Z. W., Cohen, E. G. D., R e t h m e i e r , B. C., De Boer, J., Physica g~I (1957) 407. R o s e n b l u t h , M. N. and R o s e n b l u t h , A. W., J. chem. Phys. 22 (1954) 881. M e t r o p o l i s , N., R o s e n b l u t h , A. W., R o s e n b l u t h , M. N., Teller, A. H. and Teller, E., J. chem. Phys. 21 xl (1953) I087.
Physica XXV