Prediction of a structural transition in the hard disk fluid
THE JOURNAL OF CHEMICAL PHYSICS 133, 164507 共2010兲
Prediction of a structural transition in the hard disk fluid Jarosław Piasecki,1,a兲 Piotr Szymczak,1 and John J. Kozak2 1
Institute of Theoretical Physics, University of Warsaw, Hoza 69, 00-681 Warsaw, Poland DePaul University, 243 South Wabash Avenue, Chicago, Illinois 60604-2301, USA
2
共Received 12 July 2010; accepted 30 August 2010; published online 25 October 2010; publisher error corrected 30 December 2010兲 Starting from the second equilibrium equation in the BBGKY hierarchy under the Kirkwood superposition closure, we implement a new method for studying the asymptotic decay of correlations in the hard disk fluid in the high density regime. From our analysis and complementary numerical studies, we find that exponentially damped oscillations can occur only up to a packing fraction ⴱ ⬃ 0.718, a value that is in substantial agreement with the packing fraction, ⬃ 0.723, believed to characterize the transition from the ordered solid phase to a dense fluid phase, as inferred from Mak’s Monte Carlo simulations 关Phys. Rev. E 73, 065104 共2006兲兴. Next, we show that the same method of analysis predicts that the exponential damping of oscillations in the hard sphere fluid becomes impossible when = 4n3关1 + H共1兲兴 ⱖ 34.81, where H共1兲 is the contact value of the correlation function, n is the number density, and is the sphere diameter in exact agreement with the condition, ⱖ 34.8, which is first reported in a numerical study of the Kirkwood equation by Kirkwood et al. 关J. Chem. Phys. 18, 1040 共1950兲兴. Finally, we show that our method confirms the absence of any structural transition in hard rods for the entire range of densities below close packing. © 2010 American Institute of Physics. 关doi:10.1063/1.3491039兴 I. INTRODUCTION
The second equilibrium equation in the BBGKY hierarchy establishes an exact relation between the pair and triplet number density. Invoking the Kirkwood superposition approximation yields a nonlinear integral equation for the pair correlation function.1,2 Interest in studying the analytic and numerical properties of the resulting Yvon–Born–Green and/or Kirkwood equation began with Kirkwood and co-workers,3–6 and continues to the present day.7 Of particular interest is whether the closed equation provides an essentially correct description of the fluid phase, and whether a 共possible兲 change in the analytic character of the solutions signals a change from the fluid phase to a solid phase. It is to the latter question that the methods of the present contribution are directed. One approach to explore analytically the possibility of a phase transition from the fluid phase to the solid phase is to mobilize the theory of nonlinear integral equations, focusing on theorems that establish the necessary and sufficient conditions for the existence and uniqueness of solutions, and bifurcation points.8–12 An alternative approach is to introduce a moment expansion by means of which the Yvon–Born– Green 共YBG兲 equation can be cast into a nonlinear differential equation which may be used to analyze long-range correlations.13,14 The present contribution is centered on a new method of studying the asymptotic decay of correlations, which was first introduced in Ref. 15 for the hard sphere fluid. In this paper, we focus on the hard disk fluid. The method leads to the prediction of a structural transition in a兲
Electronic mail: [email protected].
0021-9606/2010/133共16兲/164507/10/$30.00
both the hard sphere and hard disk fluids, and no transition in the hard rod system 共as must be the case兲. We shall show that the values of packing fractions at which the predicted transitions occur are in agreement with estimates derived from numerical solution of the Kirkwood equation6 and recent Monte Carlo simulations.16 II. THE SECOND EQUILIBRIUM HIERARCHY EQUATION FOR HARD DISKS
We consider a gas of hard disks of diameter at thermal equilibrium with constant number density n and temperature T. Let n2共r12兲 denote the number density of pairs of particles situated at distance r12 = 兩r1 − r2兩. Using the fact that n2共r12兲 = 0 for r12 ⬍ , we introduce a dimensionless function y 2共r12兲 defined by n2共r12兲 = n2共r12 − 兲y 2共r12兲,
共1兲
where is a unit step function. We assume that n2共r12兲 → n2 when r12 → ⬁. The twoparticle dimensionless correlation function h2共r12兲 is then defined by the cluster decomposition y 2共r12兲 = 1 + h2共r12兲, so that n2共r12兲 = 共r12 − 兲关1 + h2共r12兲兴, n2
共2兲
h2共r12兲 is thus supposed to satisfy the asymptotic condition lim h2共r12兲 = 0.
r12→⬁
共3兲
The second equilibrium YBG hierarchy equation establishes an exact relation between n2共r12兲 and the reduced three-particle number density n3共r1 , r2 , r3兲. Introducing the excluded volume factor we write the three-particle density as
133, 164507-1
© 2010 American Institute of Physics
164507-2
J. Chem. Phys. 133, 164507 共2010兲
Piasecki, Szymczak, and Kozak
n3共r1,r2,r3兲 = 共r12 − 兲共r13 − 兲 ⫻共r23 − 兲n3y 3共r12,r13,r23兲,
共4兲
The function n3 depends only on the distances rij = 兩rij兩 = 兩ri − r j兩, and is a symmetric function of the three variables r12 , r13 , r23. In the case of hard disks y 2 is related to y 3 through the second YBG hierarchy equation 共see Appendix A兲 d y 2共r兲 = n dr
冕
Our aim is to derive from Eq. 共10兲 the equivalent integral equation satisfied by the dimensionless correlation function
To this end, we insert Eq. 共11兲 into Eq. 共10兲 finding d ln关1 + H共x兲兴 dx
ˆ 共rˆ · ˆ 兲共兩r − ˆ 兩 − 兲 d
ˆ 兩兲, ⫻y 3共r, ,兩r −
= n2关1 + H共1兲兴
共5兲
ˆ and rˆ are unit vectors 兩 ˆ 兩 = 兩rˆ 兩 = 1. Putting rˆ · ˆ where = cos we rewrite Eq. 共5兲 in an explicit form d y 2共r兲 = n dr
冕
2
+
d cos y 3共r, , 冑r2 − 2r cos + 2兲 共6兲
When writing Eq. 共6兲, the equality ˆ 兩 − 兲 = 共r − 2 cos兲, 共兩r −
冕
冕
2
2
冎
⬁
dz共2 − z兲
x
冑 冉冊 冋 冑 冉冊 1−
= 共2 − x兲 −
x 2
where
x 2
2
z 2
2
1−
x 2
2
+ arccos
册
x , 2
共14兲
冋 冑 冉冊 x 2
1−
x 2
2
+ arccos
册
x , 2
共15兲
共8兲
Iⴱⴱ共x兲 = −
1 2
冕 冕 ⬁
2
dz
x
d
0
⫻cos 共z − 2 cos 兲H共冑z2 + 1 − 2z cos兲. 共16兲
共9兲
we find that it satisfies the nonlinear equation
冕
2
d cos Y共冑x2 − 2x cos + 1兲
It turns out that the angular integration in the formula for Iⴱⴱ共x兲 can be explicitly performed 共the calculation is presented in Appendix B兲. One eventually finds
0
⫻共x − 2 cos兲,
,
and
It is convenient to rewrite Eq. 共8兲 using the dimensionless distance x = r / . Denoting by Y共x兲, the function Y共x兲 = y 2共x兲,
1−
ln关1 + H共x兲兴 = 2n2关1 + H共1兲兴关Iⴱ共x兲 + Iⴱⴱ共x兲兴,
Iⴱ共x兲 = 共2 − x兲 −
⫻共r − 2cos兲.
冑 冉冊
we get
共7兲
dcos y 2共冑r2 − 2r cos + 2兲
共12兲
d共x − 2 cos兲cos = − 2共2 − x兲
0
d ln Y共x兲 = n2Y共1兲 dx
d cos H共冑x2 − 2x cos + 1兲
共13兲
d y 2共r兲 = ny 2共r兲y 2共兲 dr ⫻
2
0
III. THE KIRKWOOD SUPERPOSITION APPROXIMATION
Adopting Eq. 共7兲 leads to a closed equation
d cos 共x − 2 cos兲
0
In order to derive an integral equation for H共x兲 we integrate Eq. 共12兲 over the spatial interval 共x , ⬁兲. Using the formulae
冕
The Kirkwood superposition approximation consists in replacing in Eq. 共6兲 the three-particle density by the product of two-particle densities corresponding to three different pairs of particles
2
⫻共x − 2 cos兲 .
has been used. The rigorous relation 共6兲, which is valid for r ⬎ , will be the starting point for subsequent considerations.
ˆ 兩兲 → y 2共r兲y 2共兲y 2共兩r − ˆ 兩兲. y 3共r, ,兩r −
冕
再冕
0
0
⫻共r − 2 cos兲.
共11兲
H共x兲 = Y共x兲 − 1.
共10兲
valid in the region x ⱖ 1. Equation 共10兲 represents the closure of the YBG hierarchy corresponding to the superposition approximation.
Iⴱⴱ共x兲 = −
冕
x+1
x−1
ds共s − 1兲sH共s兲arccos
冉
冊
x2 + s2 − 1 . 2xs 共17兲
In this way, we arrive at an integral equation
164507-3
J. Chem. Phys. 133, 164507 共2010兲
Structural transition in the hard disk fluid
H
becomes slower as the surface fraction is increased, and a pronounced peak structure appears.
2
IV. LINEARIZATION OF KIRKWOOD’S EQUATION: THE RING APPROXIMATION 1
x 2
4
6
8
FIG. 1. Pair correlation function H共x兲 for the surface fraction = n2 / 4 = 0.35 共dashed line兲 and = 0.59 共solid line兲.
ln关1 + H共x兲兴 = 2n2关1 + H共1兲兴
再 冋
⫻ 共2 − x兲 − −
冕
x 2
冑 冉冊 x 2
1−
2
+ arccos
x 2
册
x+1
ds共s − 1兲sH共s兲
x−1
⫻arccos
冉
x2 + s2 − 1 2xs
冊冎
y 3共r12,r13,r23兲 = 1 + h2共r12兲 + h2共r13兲 + h2共r23兲
共19兲
where L is the integral operator given by L共H共x兲兲 ⬅ − 1 + exp 2n2关1 + H共1兲兴
再 冋
−
冕
冑 冉冊 x 1− 2
2
x + arccos 2
册
ds共s − 1兲sH共s兲
⫻arccos
冉
x2 + s2 − 1 2xs
冊冎 冊
y 3共r12,r13,r23兲 ⯝ 1 + h2共r12兲 + h2共r13兲 + h2共r23兲.
It is important here to note that while rejecting in the cluster decomposition 共22兲 the three-particle correlations we nevertheless retain in the hierarchy equation, the full excluded volume factor represented by the product of unit step functions 关see definition 共4兲兴. In fact, this factor represents the exact lowest order term in the density expansion of the threeparticle number density, and is of fundamental importance for the correct description of hard disks at low densities. Comparing Eq. 共23兲 with the superposition formula y 3共r12,r13,r23兲 ⯝ 关1 + h2共r12兲兴关1 + h2共r13兲兴关1 + h2共r23兲兴, 共24兲
.
共20兲
再
d H共x兲 = n2 − 2共2 − x兲 dx
The above integral equation for H共x兲 was solved by a standard Neumann method with successive over-relaxation.17 The iterative solutions are then given by Hn = 共1 − ␣兲Hn−1 + ␣L共Hn−1兲.
共23兲
we see that the ring approximation corresponds exactly to the linearization of the Kirkwood theory in h2. The linearized form of Eq. 共12兲 reads
x+1
x−1
共22兲
Neglecting h3 is thus equivalent to the approximation
H共x兲 = L共H共x兲兲,
x ⫻ 共2 − x兲 − 2
+ h3共r12,r13,r23兲.
共18兲
,
representing the superposition closure for the correlation function H共x兲. Our method of determining H共x兲 is based on the fact that the solution of Eq. 共18兲 satisfying the boundary condition limx→⬁ H共x兲 = 0 can be obtained numerically by iterations. Note that we can rewrite Eq. 共18兲 in the following form:
冉
Before continuing the analysis based on Eq. 共18兲, let us make a comment on the relationship between the superposition approximation and the ring approximation, well known from the kinetic theory 共see Ref. 15 and references given therein兲. Originally, the ring approximation was applied to the study of long wavelength hydrodynamic phenomena, and was defined by neglecting in the second equation of the dynamical BBGKY hierarchy all contributions from the threeparticle correlations. The three-particle correlation function h3 is defined by the cluster decomposition of the density y 3
+
2
1−
x 2
2
关1 + H共1兲 + H共x兲兴
dcos 共x − 2 cos兲
0
冎
⫻H共冑x2 + 1 − 2x cos兲 .
共21兲
The relaxation parameter ␣ was taken to be 0.25. The iterations were continued until successive values of H共x = 1兲 differed by less than ⑀ = 10−5, except in the vicinity of the threshold surface fraction ⴱ 共see Sec. V B兲, where the convergence was slow and the iterations were discontinued at ⑀ = 10−2. Examples of the correlation functions obtained in this way are presented in Fig. 1. As it is seen, the decay of H共x兲
冕
冑 冉冊
共25兲
In order to derive an integral equation for H共x兲, we proceed as before by integrating Eq. 共25兲 over the spatial interval 共x , ⬁兲. The result reads H共x兲 = 2n2关Jⴱ共x兲 + Jⴱⴱ共x兲兴, where
164507-4
J. Chem. Phys. 133, 164507 共2010兲
Piasecki, Szymczak, and Kozak
Z
Z共兲 =
10
p = 1 + n2共1 + H共1兲兲 = 1 + 2共1 + H共1兲兲, nkT 2 共30兲
where the contact value H共1兲 is also a function of . For comparison, we include here Z共兲 dependence as predicted by the scaled particle theory 共SPT兲18
5
0
0.1
ZSPT共兲 = x
0 0.2
0.3
0.4
0.5
0.6
FIG. 2. Compressibility factor Z vs the surface fraction for the full Kirkwood approximation 共circles兲, the ring approximation 共squares兲, and the scaled particle theory 关Eq. 共31兲, dashed line兴.
Jⴱ共x兲 = 共2 − x兲
再冋
冑 冉冊 册 冕 冑 冉冊 冎
x 2
−
2
1−
x 2
dz
1−
2
⫻关1 + H共1兲兴 +
x
+ arccos
x 2
关H共z兲兴 ,
共26兲
and Jⴱⴱ共x兲 = −
1 2
冕 冕 ⬁
2
dz
d cos 共z − 2 cos兲
0
x
⫻H共冑z2 + 1 − 2zcos兲 = Iⴱⴱ共x兲, 关see Eq. 共16兲兴. Using again the calculation presented in Appendix B, we eventually find Jⴱⴱ共x兲 = −
冕
x+1
ds共s − 1兲sH共s兲arccos
x−1
冉
冊
x2 + s2 − 1 . 2xs 共27兲
H共x兲 = 2n2共2 − x兲
再冋
冑 冉冊 冕 冑 冉冊 x 1− 2
2
⫻关1 + H共1兲兴 + − 2n2
冕
dz
x
1−
z 2
2
x + arccos 2
2
关H共z兲兴
x+1
x−1
ds共s − 1兲sH共s兲arccos
As it is seen, the iteration results agree with the SPT predictions up to approximately = 0.4.19 For larger packing fractions, the Kirkwood approximation tends to underestimate the compressibility factor, whereas the ring approximation overestimates it.
A. Breakdown of the method used for attractive interactions
The fundamental information concerning the internal structure of the system is contained in the spatial dependence of correlations. In particular, the law governing the asymptotic vanishing of correlations is of primary importance. In order to determine the behavior of H共x兲 at large distances, it seems natural to follow the moment analysis presented in Refs. 13 and 14. The calculation would proceed as follows. For x ⬎ 2, the Kirkwod equation 共12兲 can be conveniently written as d ln关1 + H共x兲兴 = n2关1 + H共1兲兴 dx
ˆ 共xˆ · ˆ 兲H共兩x − ˆ 兩兲. d
冉
2
冎
册
ˆ 兩兲 around When x Ⰷ 1, the power series expansion of H共兩x − the point 兩x兩 = x yields nonzero contributions only from terms ˆ 兲. The calculation up to involving odd powers of cos = 共xˆ · the third derivative of H yields the expansion ln关1 + H共x兲兴⬘ = n2关1 + H共1兲兴
2
冊
⫻
x +s −1 , 2xs
representing the ring approximation. Equation 共28兲 can again be solved numerically by iterations, and is expected to yield relevant results in the low density regime. The comparison of the results obtained with the full Kirkwood approximation and its linearization is presented in Fig. 2. Defining the surface fraction as
2 , 4
we plot here the compressibility factor Z共兲 given by
共29兲
冕
2
再
cos − cos H⬘共x兲
0
冉 冊
1 H⬘共x兲 ⬘ − 关cos − 共cos兲3兴 2 x
共28兲
= n
冕
共32兲
In this way we arrive at an integral equation x − 2
共31兲
V. ASYMPTOTIC DECAY OF CORRELATIONS: PREDICTING A STRUCTURAL TRANSITION
2
z 2
1 . 共1 − 兲2
冎
1 − 共cos兲3H共x兲 + ¯ , 6
共33兲
where ⬘ denotes the derivative with respect to x. Using the boundary condition limx→⬁ H共x兲 = 0 together with the equalities
冕
2
0
d共cos兲2 = ,
冕
2
0
3 d共cos兲4 = , 4
and adopting for large distances the asymptotic formula
164507-5
J. Chem. Phys. 133, 164507 共2010兲
Structural transition in the hard disk fluid
ln关1 + H共x兲兴 ⬃ H共x兲,
ReΚ,ImΚ 5
we find a linear differential equation of the form 1 H⬙共x兲 + H⬘共x兲 + ␣2H = 0, x
共34兲
A
0 4
where
冉
␣2 = 8 1 +
冊
1 ⬎ 0. n2关1 + H共1兲兴
B. Prediction of a structural transition
Let us consider again the region x ⬎ 2 where the Eq. 共32兲 holds. As limx→⬁ H共x兲 = 0, we can replace in Eq. 共32兲 the function ln关1 + H共x兲兴 by H共x兲, and consider the equation
冕
2
d cos H共冑x2 − 2x cos + 1兲,
12
16
20
5
However, Eq. 共34兲 is the equation for the Bessel function J0共␣x兲. The solution of Eq. 共34兲 is thus an oscillating function whose amplitude decays as 1 / 冑x. Unfortunately, the above result is in disagreement with numerical predictions of exponentially damped oscillations. Use of the moment expansion as developed in Ref. 13 leads to erroneous results. Clearly one has to take into account the ˆ 兩兲 to get whole infinite series in the expansion of H共兩x − reliable predictions. We have thus to give up this kind of expansion, and look for a different approach. In a broad outline, the failure of the moment expansion developed in Refs. 13 and 14 to describe the structural transition in the hard disk fluid can be traced to the nature of the governing intermolecular potential. The studies13,14 deal with the description of correlations in the vicinity of the liquidvapor critical point, where attractive forces play a dominant role. Using the moment expansion, one recovers at large distances the classical Ornstein–Zernike formula. The present study deals with the fusion transition, where short-range, repulsive forces play the critical role. The new method introduced here 共see Sec. V B兲 effectively accounts for the difference in the potential governing these two transitions and, anticipating our later results, leads to an analytic prediction of the packing fraction at which a structural transition occurs in the hard disk 共and hard sphere兲 fluid.
A d H共x兲 = dx 2
8
共35兲
0
FIG. 3. The real part 共solid line兲 and imaginary part 共dashed line兲 of the root of Eq. 共38兲 corresponding to the slowest decaying mode.
notice that all derivatives of H共x兲 satisfy the same equation. It is thus natural to consider H共x兲 as a linear combination of exponential modes exp共x兲, where is a complex number. The function exp共x兲 satisfies Eq. 共37兲 provided solves the equation
=
A 2
冕
2
d cos exp共− cos兲 = − AI1共兲,
共38兲
0
where I1共兲 is a modified Bessel function. The physically acceptable solutions = a + ib are those with a negative real part a ⬍ 0 that assures the exponential damping of oscillations. The first root 共with the smallest absolute value of a兲 corresponding to the slowest decay of the correlation function is shown in Fig. 3. The values of 共兲 predicted with the use of Eq. 共38兲 are in good agreement with the decay of the amplitude of H共x兲 determined by the iterative solution of integral equation 共19兲. Plotting the absolute value of H共x兲 on a logarithmic plot, and fitting it by the single mode ␣e⬘x, we obtain values of ⬘ that are slightly below those obtained by solving Eq. 共38兲. For example, for = 0.3 we get ⬘ ⬇ −2.1 共cf. Fig. 4兲, whereas the corresponding value of for that surface fraction is ⬇ −2.02. Similarly, for = 0.55, we get, respectively, ⬘ ⬇ −0.8 and ⬇ −0.7. The fact that the values of ⬘ remain slightly below those of can be understood by noting that corresponds to the slowest decaying mode, whereas in the numerical data on H共x兲 we also see nonzero contributions from other, faster decaying modes. An interesting feature of the 共兲 dependence presented in Fig. 3 is the fact that the root becomes purely imaginary at
where H 1
A = 2n2关1 + H共1兲兴. We then use the expansion 2 冑x2 − 2x cos + 1 = x − cos + sin + ¯ , 2x
0.01
共36兲
106
to arrive at the equation A d H共x兲 = dx 2
冕
2
d cos H共x − cos兲,
104
108
共37兲
1010
0
valid for x Ⰷ 1. In order to determine the large x behavior of correlations, we have thus to analyze the solution of Eq. 共37兲. We
2
4
6
8
10
x
FIG. 4. The absolute value of H共x兲 for = 0.3 together with a corresponding fit of the form ae⬘x with ⬘ = −2.1.
164507-6
J. Chem. Phys. 133, 164507 共2010兲
Piasecki, Szymczak, and Kozak
8 ⬎ 15.12,
A 16
13
12
共41兲
representing for a physically impossible condition 共beyond close packing兲. The ring approximation is unable to describe a qualitative change in the hard disk correlation function. For any accessible surface fraction, it predicts exponentially damped oscillations.
15 14
⬎ 1.89,
or
11
10
C. Structural transition in a hard sphere fluid
x 0.5
0.52
0.54
0.56
0.58
0.6
0.62
FIG. 5. The extrapolation of A共兲 dependence. The points correspond to the values of A obtained from the iterative procedure, and the dashed line is given by A = Aⴱ ⬇ 15.12.
A = Aⴱ ⬇ 15.1. This point can be made more precise by analyzing when Eq. 共38兲 acquires a purely imaginary solution = ib. As I1共ib兲 = iJ1共b兲, Eq. 共38兲 implies the condition 1 J1共b兲 1 = 关J0共b兲 + J2共b兲兴 = − , 2 b A
共39兲
where the first equality follows from the recurrence relation between Bessel functions Jn. The absolute minimum of the sum of Bessel functions 关J0共b兲 + J2共b兲兴 equals ⫺0.1323. The necessary condition for the disappearance of damping has thus the form 2 ⱕ 0.1323, A
or
2 = 4n3
冉
共40兲
The above inequality shows that the nature of correlations could change for sufficiently high values of the surface fraction 关Eq. 共29兲兴 occupied by hard disks. The determination of the value of ⴱ corresponding to the equality Aⴱ = 8ⴱ关1 + H共x = 1; ⴱ兲兴 = 15.12, requires the knowledge of the contact value H共1兲 as a function of the surface fraction. We have studied this question numerically. An accurate estimation of the precise value of ⴱ corresponding to Aⴱ is difficult because as we approach ⴱ, the iteration procedure demands an increasingly larger number of iterations to converge to a solution with the required accuracy. Additionally, a computational domain over which the solution is sought must also be progressively extended as we approach ⴱ, since the decay of H共x兲 is very weak there. We estimated the value of ⴱ by calculating A共兲 for several values of in the range of 0.5ⱕ ⱕ 0.6 and then extrapolating to larger values of , as illustrated in Fig. 5. In this way, we obtain the estimate of ⴱ ⬇ 0.622. Let us close this section with a comment on the ring approximation. In order to find the asymptotic behavior of correlations satisfying Eq. 共28兲, it is sufficient to consider the linearized version of Eq. 共35兲. However, this is equivalent to the replacement of the factor A = 2n2关1 + H共1兲兴 by 2n2 = 8. The inequality 共40兲 becomes then
冊
sh共兲 − ch共兲 .
共42兲
In passing to the Kirkwood superposition approximation, we need only to replace the factor 4n3 in the above equation by = 4n3关1 + H共1兲兴. This fact follows from the previously made remark that the ring approximation represents exactly the linearized Kirkwood theory. In order to determine the range of values of ⬎ 0 beyond which the exponential damping becomes impossible we explore the possibility of vanishing of the real part a of the complex number = a + ib. Equation 共42兲 reduces then to
冉
b2 = cosb −
A = 2n2关1 + H共1兲兴 ⱖ Aⴱ
Aⴱ = 15.1171 ¯ ⯝ 15.12.
The integral equation satisfied by the correlation function H共x兲 of a three-dimensional hard sphere fluid for x Ⰷ 1 has been analyzed in Ref. 15 within the ring approximation. The frequencies of the exponential modes exp共x兲 describing the long distance decay of correlations are in this case solutions of the equation 关see Eq. 共22兲 in Ref. 15兴
冊
sin b , b
or
=
b3 . b cosb − sin b
共43兲
The absolute minimum of the function y共b兲 =
b3 , b cosb − sin b
共44兲
in the region where y共b兲 ⬎ 0 equals 共45兲
y min = 34.81 ¯ .
Hence, for ⬎ 34.81, Eq. 共42兲 acquires purely imaginary solutions and the exponential damping vanishes. It is quite remarkable that 60 years ago, Kirkwood et al.6 concluded from numerical studies of the integral equation for H共x兲 that when exceeds 34.8 the correlation function H共x兲 is not integrable anymore. Our method provides a simple analytic confirmation of this result.
D. The hard rod fluid: Testing the method
The rigorous calculation of the two-particle correlation function for hard rods 共see e.g., Ref. 20兲 shows that its structure corresponds to exponentially damped oscillations at all possible densities. One finds ⬁
n关H共x兲 + 1兴 = 兺 关x − 共k + 1兲兴 k=0
k关x − 共k + 1兲兴k k!
⫻exp兵− 关x − 共k + 1兲兴其, where
共46兲
164507-7
J. Chem. Phys. 133, 164507 共2010兲
Structural transition in the hard disk fluid
Κ ns
0 0.2
0.4
0.6
0.8
1
H 5
2
4
4
3 2
6
1
x
8
2
FIG. 6. The real part of n = − + Wn共 exp 兲, n = 1 , . . . , 5 共top to bottom兲 as a function of n = / 共1 + 兲.
=
n . 1 − n
In fact, the superposition law turns out to be exact for a one-dimensional hard rod fluid,20 and the second equation of the equilibrium hierarchy takes a particularly simple form H⬘共x兲 = n关共x − 2兲H共x − 1兲 − H共x兲 − 共2 − x兲兴 ⫻关H共1兲 + 1兴.
共47兲
The formula 共46兲 represents the solution of Eq. 共47兲. When x ⬎ 2, we find H⬘共x兲 = n关H共x − 1兲 − H共x兲兴关H共1兲 + 1兴.
共49兲
It turns out that all solutions of Eq. 共49兲 can be expressed in terms of the multivalued Lambert W function. Indeed, the special function W共z兲 is defined on the complex plane by the equation z = W共z兲exp关W共z兲兴.
共50兲
But, Eq. 共49兲 can be rewritten as 共 + 兲exp共 + 兲 = exp .
共51兲
It follows that
= − + W共 exp 兲.
6
8
0 FIG. 7. The exact form of the correlation function for hard-rod fluid 共solid兲 at n = 0.8 and its asymptotic form H共x兲 = Aex with 共兲 calculated using Eq. 共52兲.
As illustrated in Figs. 7 and 8, the correlation function H共x兲 rapidly approaches the asymptotic form given by the slowest decaying mode Ae1x. Note that not only the exponential decay rate but also the oscillation period of the function agree with that given by Ae1x starting from x ⬇ 3. This shows that the other modes play a negligible role in influencing the behavior of H共x兲 for intermediate and large x values, thus lending further support to our approach of focusing on the slowest decay mode only.
共48兲
Applying the same method as that used for hard disks, we look for exponential modes exp共x兲 solving Eq. 共48兲. The complex frequency satisfies the equation
= n关H共1兲 + 1兴共e− − 1兲 = 共e− − 1兲.
4
共52兲
The principal branch of Lambert function W0共z兲 obeys W0共xex兲 = x for real x, thus = 0 for that branch. However, other branches Wn共z兲 with n ⬎ 0 give values of n with a negative real part corresponding to the decay of the correlation function. The first five solutions, n共兲 , n = 1 , . . . , 5, are presented in Fig. 6 as functions of n = / 共1 + 兲 共modes with larger n decay faster兲. The negative real part of n , n ⬎ 0 is different from zero for any accessible density, and vanishes only at close packing where n = 1. There are thus no purely imaginary solutions = ib of Eq. 共49兲. This confirms the correctness of the method, showing the absence of any structural transition in one dimension over the whole range of densities below close packing.
VI. DISCUSSION AND CONCLUSIONS
The question of whether a system of particles interacting via a purely repulsive potential 共only兲 can undergo a phase transition has been under continuous investigation since first posed and addressed by Kirkwood over 70 years ago. For a system of hard disks, the first numerical evidence was provided by Alder and co-workers.21,22 The data reported in Ref. 21 showed that the hard disk freezing transition occurred at a density smaller than the density of closest packing 关corresponding to an area fraction of 0 = / 冑12= 0.906 90兴, and suggested that the liquid to solid transition was first order. The most recent Monte Carlo simulations 共on a system of 4 ⫻ 106 disks兲 of Mak16 suggest that melting consists of a continuous transition from the ordered solid to an intermediate 共hexatic兲 phase23–27 at a packing fraction = 0.723, and
H 1 0.1 0.01 0.001 104 105 2
4
6
8
10
12
14
x
FIG. 8. The absolute value of H共x兲 − 1 for hard rod fluid 共solid兲 at n = 0.8 and the exponential asymptote H共x兲 = AeRe共兲x.
164507-8
J. Chem. Phys. 133, 164507 共2010兲
Piasecki, Szymczak, and Kozak
either a very weak first-order or a continuous transition from the intermediate phase to the fluid phase at a packing fraction = 0.699. From the analysis and numerical evidence presented in Sec. V B, we have calculated the area fraction at which a structural change in the hard disk fluid can take place, viz., ⴱ ⬃ 0.622. Converting this area fraction to a packing fraction gives ⴱ = ⴱ / 0 ⬃ 0.718. When compared to the estimate reported by Mak for the transition from the ordered solid phase to a dense fluid phase, ⬃ 0.723, one finds the two values are in substantial agreement. Also of interest is our prediction of a structural transition in a system of hard spheres. As noted in Sec. V C, Kirkwood, Maum and Alder6 found that for values of ⱖ 34.8, no solutions of the YBG and Kirkwood integral equations exist for which, in their notation, x2关g共x兲 − 1兴 is integrable. We find that the value of beyond which exponential damping becomes impossible is = 34.81. Hence, the results for hard spheres appear to be in exact agreement. When we consider a system of hard rods, the analytic method developed here confirms the absence of any structural transition in d = 1 for the entire range of densities below close packing 关a result that was already known to Rayleigh 共Ref. 28兲 and Korteweg 共Ref. 29兲兴. Our prediction of a structural phase transition is based on the analysis of an integral equation whose derivation assumed sufficiently fast decay of correlations. If a phase becomes ordered and correlations do not decay, the integral equation to which our method was applied does not hold. This can preclude the possibility of studying a region where a new equilibrium phase may be formed. In particular, the identification of a region intermediate between melting and freezing 共see following paragraph兲 and the characterization of an ordered 共solid兲 phase is certainly beyond the scope of our approach. To study these questions within the Kirkwood superposition approximation, one must go back to the original YBG hierarchy equation 共Kirkwood’s closure does not assume the rapid decay of correlations兲. To develop this point further, there are three structural aspects of the hard disk transition that are not captured by the method developed in this paper. First, we find no evidence for the existence of an intermediate, or hexatic, phase predicted by the Kosterlitz, Thouless, Halperin, Nelson, and Young 共KTHNY兲 theory of d = 2 melting,23–27 and supported by Mak’s Monte Carlo simulations. Second, we find no evidence for the development of a shoulder on the second maximum of disk radial distribution function in the vicinity of the freezing transition 共 = 0.686兲 reported by Truskett et al.30 based on molecular dynamics simulations, and later correlated with structural rearrangements occurring at increasing disk density.31–34 Third, we cannot confirm the existence of regions of fivefold coordination in the dense fluid phase, first predicted by Bernal35–37 based on his “ball and spoke” model of a random assembly of hard-core particles, although it has been conjectured that the hexatic phase might be correlated with randomly dispersed regions of 5-, 6-, and 7-member disk clusters forming percolated tessellations that span the d = 2 space.38 The larger point, however, relates to the original Kirkwood prediction, viz., that a system of particles interacting
via purely repulsive forces, here hard disks but also hard spheres, can undergo a phase transition. Although not widely accepted at first, following the work of Onsager on the isotropic-nematic transition in a d = 3 dimensional system of thin hard rods,39 there developed a gradual realization that a phase transition can be entropy driven. As elaborated by Frenkel,40 in hard-core systems the entropy in the ordered crystalline phase is actually larger than in the fluid phase; quoting directly, “the entropy decreases because the density is no longer uniform in orientation or position, but the entropy increases because the free-volume per particle is larger in the ordered than in the disordered phase.” The present contribution provides further evidence for the essential correctness of Kirkwood’s insight.
APPENDIX A: DERIVATION OF EQ. „5… FROM THE BBGKY HIERARCHY
Consider a gas of hard disks of mass m and diameter . We denote by j ⬅ 共r j , v j兲, j = 1 , 2 , . . ., the one-particle state in which a disk has position r j and velocity v j. The average number density of s-particle clusters occupying at time t the s-particle state 共1 , 2 , . . . , s兲 is called the s-particle reduced distribution f s共1 , 2 , . . . , s ; t兲. The dynamical evolution of the hard disk fluid is described in the thermodynamic limit by the BBGKY hierarchy equations. The second of them establishes a relation between f 2 and f 3
冉
冊
+ v1 · + v2 · − ¯T共1,2兲 f 2共1,2;t兲 t r1 r2 =
冕
¯ 共1,3兲 + ¯T共2,3兲兲f 共1,2,3;t兲. d3共T 3
共A1兲
The effects of binary collisions are described by the operator ¯T共i , j兲 ¯T共i, j兲 =
冕
ˆ v12 · ˆ 共vij · ˆ 兲关␦共rij − 兲bˆ d
− ␦共rij + 兲兴.
共A2兲
Here the Dirac ␦-distributions restrict the distances 兩rij兩 between the centers of the disks at the moment of impact to ˆ is oriented perpentheir diameter = 兩兩. The vector = dicularly to the surface of colliding disks at the point of impact. The action of the operator bˆ consists in replacing the velocities vi , v j by their precollisional values v⬘i , v⬘j corresponding to the inverse elastic collision. As in elastic collisions the kinetic energy is conserved, in the case of products of Maxwell distributions
共v兲 =
冉
冊
m exp共− mv2/2kBT兲, 2 k BT
we find bˆ 关共vi兲共v j兲兴 = 共v⬘i 兲共v⬘j 兲 = 共vi兲共v j兲. Hence, in the case of equilibrium reduced distributions
共A3兲
164507-9
J. Chem. Phys. 133, 164507 共2010兲
Structural transition in the hard disk fluid
f s共1,2, . . . ,s兲 = ns共r1,r2 ¯ rs兲共v1兲共v2兲 ¯ 共vs兲, 共A4兲 Eq. 共A1兲 takes the form
冉
v1 · + v2 · − r1 r2
冕
ˆ v12 · ˆ ␦共rij − 兲 d
冊
冕
dr3dv3
冕
− r12
冕
冕 冕 dr3
共A5兲
冊
Iⴱⴱ共x兲 = −
冕
ˆ v12 · ˆ ␦共r12 − 兲, d
y 2共r12兲. r12
冕
x−cos
ds共s − cos兲 共B2兲
d
ds共x + 冑1 − 2 − s兲共s − x + 冑1 − 2兲
⫻H共冑s2 + 2兲共s − 冑1 − 2兲.
共A7兲
共B3兲
We now introduce a new integration variable
共A6兲
z = 冑s2 + 2 . As zdz = sds, we find
冕 冕 1
H共z兲共冑z2 − 2 2 冑 z − 0 − 冑1 − 2兲共冑1 − 2 − 兩x − 冑z2 − 2兩兲
Iⴱⴱ共x兲 = − 共A8兲
=−
d
冕 冕
dz
z
2
dz zH共z兲共z − 1兲 1
d
⫻
1
冑z2 − 2 共
冑1 − 2 − 兩x − 冑z2 − 2兩兲.
ˆ v12 · ˆ 共兩共r12 − 兲兩 − 兲 d
共B4兲 共A9兲
As the equality L = R must hold for any value of the relative velocity v12, we finally find 共when r12 ⱖ 兲
冕
冕 冕 1
0
⫻y 3共r12, ,兩共r12 − 兲兩兲.
dy 2共r12兲 = n dr12
x+cos
0
On the right-hand side owing to the presence of ␦-distributions, we can perform integration over variable r3, thus obtaining R = n共r12 − 兲
冕
d cos
Putting then = sin , we get
which reduces the left hand side of Eq. 共A6兲 to L = 共r12 − 兲v12 ·
/2
⫻H共冑s2 + sin2 兲.
We now use the identity
共r12 − 兲 ⬅ r12
冕
0
⫻共r23 − 兲 − 兲y 3共r12, ,r23兲
v12 ·
dcos 关共z − 2 cos兲
0
Changing the order of integrations with the use of the asymptotic decay of the correlation function we arrive at a convenient formula
Iⴱⴱ共x兲 = −
ˆ 关共v1 · ˆ ␦共r13 − 兲 d
ˆ ␦共r23 − 兲共r13 − 兲y 3共r12,r13, 兲兴. + v2 ·
/2
dz
共B1兲
ˆ v12 · ˆ ␦共r12 − 兲 共r12 − 兲y 2共r12兲 d
= n共r12 − 兲
冕 冕
⫻H共冑z2 + 1 − 2z cos兲 − H共冑z2 + 1 + 2z cos兲兴.
ˆ 关v13 · ˆ ␦共r13 − 兲 d
We can divide both sides of Eq. 共A5兲 by the product 共v1兲共v2兲, and perform the integration over the v3 variable. Moreover, we introduce explicitly the excluded volume factors by using Eqs. 共1兲 and 共4兲. The hierarchy equation becomes v12 ·
d cos 共z − 2 cos兲
0
x
x
⫻共v1兲共v2兲共v3兲.
2
dz
⬁
=−
ˆ ␦共r23 − 兲兴n3共r12,r13,r23兲 + v23 ·
冉
冕 冕 ⬁
⫻H共冑z2 + 1 − 2z cos兲
⫻n2共r12兲共v1兲共v2兲 =
1 2
Iⴱⴱ共x兲 = −
ˆ rˆ 12 · ˆ 共兩r12 − 兩 − 兲 d
⫻y 3共r12, ,兩r12 − 兩兲,
Here the inequality 冑1 − 2 ⬎ 兩x − 冑z2 − 2兩 is equivalent to
冑z2 − 2 ⬎ x
2
+ z2 − 1 , 2x
which leads to the formula 共A10兲
which is Eq. 共5兲 of Sec. II.
Iⴱⴱ共x兲 = −
冕 冕
dzzH共z兲共z − 1兲 1
⫻ APPENDIX B: DERIVATION OF EQ. „17…
We perform here the angular integration in the contribution to the correlation function
0
d
1
冑z2 −
2
冉冑
z2 − 2 −
冊
x2 + z2 − 1 . 2x 共B5兲
Putting = z, we find
164507-10
J. Chem. Phys. 133, 164507 共2010兲
Piasecki, Szymczak, and Kozak
Iⴱⴱ共x兲 = −
冕 冕
1/z
d
⫻
0
1
冑1 −
2
9
冉冑
1 − 2 −
冊
x2 + z2 − 1 . 共B6兲 2xz
The change of the integration variable w = 冑1 − 2 yields the formula Iⴱⴱ共x兲 = −
冕 冕 冕 冕 冕
⫻ =−
=−
⫻
dz zH共z兲共z − 1兲 1
dw
冑1−1/z2 冑1 − w
⫻
冉
w− 2
x2 + z2 − 1 2xz
1
dw
冑1 − w
冉
1− 2
x2 + z2 − 1 2xz
冉
x2 + z2 − 1 − arcsin 2 2xz
冕
冊
dz zH共z兲共z − 1兲关1 − 共x − z兲2兴
冊冎
x+1
x−1
共B7兲
.
Finally, we arrive at the result Iⴱⴱ共x兲 = −
冊
dz zH共z兲共z − 1兲
共x2+z2−1兲/2xz
再
dz zY共z兲共z − 1兲arccos
冉
冊
x2 + z2 − 1 , 2xz 共B8兲
used in Eq. 共17兲. T. L. Hill, Statistical Mechanics 共McGraw-Hill, New York, 1956兲. S.A. Rice and P. Gray, The Statistical Mechanics of Simple Liquids 共Interscience, New York, 1965兲. 3 J. G. Kirkwood, J. Chem. Phys. 3, 300 共1935兲. 4 J. G. Kirkwood and E. Monroe, J. Chem. Phys. 9, 514 共1941兲. 5 J. G. Kirkwood and E. M. Boggs, J. Chem. Phys. 10, 394 共1942兲. 6 J. G. Kirkwood, E. K. Maum, and B. J. Alder, J. Chem. Phys. 18, 1040 共1950兲. 7 J. P. Hansen and I. R. McDonald, Theory of Simple Liquids, 3rd ed. 共Academic, London, 2006兲. 1 2
J. D. Weeks, S. A. Rice, and J. J. Kozak, J. Chem. Phys. 52, 2416 共1970兲. J. J. Kozak, Adv. Chem. Phys. 40, 229 共1979兲. 10 B. Bagchi, C. Cerjan, and S. A. Rice, Phys. Rev. B 28, 6411 共1983兲. 11 A. D. J. Haymet, Annu. Rev. Phys. Chem. 38, 89 共1987兲. 12 M. Baus and C. F. Tejero, Equilibrium Statistical Mechanics: Phases of Matter and Phase Transitions 共Springer, New York, 2008兲. 13 G. L. Jones, J. J. Kozak, E. Lee, S. Fishman, and M. E. Fisher, Phys. Rev. Lett. 46, 795 共1981兲; M. E. Fisher and S. Fishman, ibid. 47, 421 共1981兲. 14 M. E. Fisher and S. Fishman, J. Chem. Phys. 78, 4227 共1983兲 共This reference provides an extensive account of the primarily analytic studies of the YBG equation using the moment expansion.兲; G. L. Jones, E. K. Lee, and J. J. Kozak, ibid. 79, 459 共1983兲 共This reference provides a similar account of the primarily numerical studies of the YBG equation.兲. 15 J. Piasecki and R. Soto, Physica A 379, 409 共2007兲. 16 C. H. Mak, Phys. Rev. E 73, 065104 共2006兲. 17 R. R. Kleinman and P. M. van den Berg, PIER 5, 67 共1991兲. 18 H. Reiss, H. L. Frisch, and J. L. Lebowitz, J. Chem. Phys. 31, 369 共1959兲. 19 K. Helfand, H. L. Frisch, and J. L. Lebowitz, J. Chem. Phys. 34, 1037 共1961兲. 20 Z. Salsburg, R. Zwanzig, and J. G. Kirkwood, J. Chem. Phys. 21, 1098 共1953兲. 21 B. J. Alder and T. E. Wainwright, Phys. Rev. 127, 359 共1962兲. 22 B. J. Alder, W. G. Hoover, and T. E. Wainwright, Phys. Rev. Lett. 11, 241 共1963兲. 23 J. M. Kosterlitz and D. J. Thouless, J. Phys. 共Paris兲 C5, 1124 共1972兲. 24 B. I. Halperin and D. R. Nelson, Phys. Rev. Lett. 41, 121 共1978兲. 25 D. R. Nelson and B. I. Halperin, Phys. Rev. B 19, 2457 共1979兲. 26 A. P. Young, Phys. Rev. B 19, 1855 共1979兲. 27 For a comprehensive review, see: D. R. Nelson, in Phase Transitions and Critical Phenomena, edited by C. Domb and J. L. Lebowitz 共Academic, London, 1983兲, Vol. 7. 28 Lord Rayleigh, Nature 共London兲 45, 80 共1891兲. 29 D. T. Korteweg, Nature 共London兲 45, 152 共1891兲. 30 T. M. Truskett, S. Torquato, S. Sastry, P. G. Debenedetti, and F. H. Stillinger, Phys. Rev. E 58, 3083 共1998兲. 31 A. Huerta, G. G. Naumis, D. T. Wasan, D. Henderson, and A. Trokhymchuk, J. Chem. Phys. 120, 1506 共2004兲. 32 F. Moucka and I. Nezbeda, Phys. Rev. Lett. 94, 040601 共2005兲. 33 A. Huerta and G. G. Naumis, Phys. Rev. Lett. 90, 145701 共2003兲. 34 A. Huerta, D. Henderson, and A. Trokhymchuk, Phys. Rev. E 74, 061106 共2006兲. 35 J. D. Bernal, Trans. Faraday Soc. 33, 27 共1937兲. 36 J. D. Bernal, Nature 共London兲 183, 141 共1959兲. 37 J. D. Bernal and S. V. King, Physics of Simple Fluids 共John Wiley & Sons, Inc., New York, 1968兲, p. 231. 38 J. J. Kozak, J. Brzezinski, and S. A. Rice, J. Phys. Chem. B 112, 16059 共2008兲. 39 L. Onsager, Ann. N.Y. Acad. Sci. 51, 627 共1949兲. 40 D. Frenkel, Physica A 263, 26 共1999兲. 8
dz zH共z兲共z − 1兲
Prediction of a structural transition in the hard disk fluid Jarosław Piasecki,1,a兲 Piotr Szymczak,1 and John J. Kozak2 1
Institute of Theoretical Physics, University of Warsaw, Hoza 69, 00-681 Warsaw, Poland DePaul University, 243 South Wabash Avenue, Chicago, Illinois 60604-2301, USA
2
共Received 12 July 2010; accepted 30 August 2010; published online 25 October 2010; publisher error corrected 30 December 2010兲 Starting from the second equilibrium equation in the BBGKY hierarchy under the Kirkwood superposition closure, we implement a new method for studying the asymptotic decay of correlations in the hard disk fluid in the high density regime. From our analysis and complementary numerical studies, we find that exponentially damped oscillations can occur only up to a packing fraction ⴱ ⬃ 0.718, a value that is in substantial agreement with the packing fraction, ⬃ 0.723, believed to characterize the transition from the ordered solid phase to a dense fluid phase, as inferred from Mak’s Monte Carlo simulations 关Phys. Rev. E 73, 065104 共2006兲兴. Next, we show that the same method of analysis predicts that the exponential damping of oscillations in the hard sphere fluid becomes impossible when = 4n3关1 + H共1兲兴 ⱖ 34.81, where H共1兲 is the contact value of the correlation function, n is the number density, and is the sphere diameter in exact agreement with the condition, ⱖ 34.8, which is first reported in a numerical study of the Kirkwood equation by Kirkwood et al. 关J. Chem. Phys. 18, 1040 共1950兲兴. Finally, we show that our method confirms the absence of any structural transition in hard rods for the entire range of densities below close packing. © 2010 American Institute of Physics. 关doi:10.1063/1.3491039兴 I. INTRODUCTION
The second equilibrium equation in the BBGKY hierarchy establishes an exact relation between the pair and triplet number density. Invoking the Kirkwood superposition approximation yields a nonlinear integral equation for the pair correlation function.1,2 Interest in studying the analytic and numerical properties of the resulting Yvon–Born–Green and/or Kirkwood equation began with Kirkwood and co-workers,3–6 and continues to the present day.7 Of particular interest is whether the closed equation provides an essentially correct description of the fluid phase, and whether a 共possible兲 change in the analytic character of the solutions signals a change from the fluid phase to a solid phase. It is to the latter question that the methods of the present contribution are directed. One approach to explore analytically the possibility of a phase transition from the fluid phase to the solid phase is to mobilize the theory of nonlinear integral equations, focusing on theorems that establish the necessary and sufficient conditions for the existence and uniqueness of solutions, and bifurcation points.8–12 An alternative approach is to introduce a moment expansion by means of which the Yvon–Born– Green 共YBG兲 equation can be cast into a nonlinear differential equation which may be used to analyze long-range correlations.13,14 The present contribution is centered on a new method of studying the asymptotic decay of correlations, which was first introduced in Ref. 15 for the hard sphere fluid. In this paper, we focus on the hard disk fluid. The method leads to the prediction of a structural transition in a兲
Electronic mail: [email protected].
0021-9606/2010/133共16兲/164507/10/$30.00
both the hard sphere and hard disk fluids, and no transition in the hard rod system 共as must be the case兲. We shall show that the values of packing fractions at which the predicted transitions occur are in agreement with estimates derived from numerical solution of the Kirkwood equation6 and recent Monte Carlo simulations.16 II. THE SECOND EQUILIBRIUM HIERARCHY EQUATION FOR HARD DISKS
We consider a gas of hard disks of diameter at thermal equilibrium with constant number density n and temperature T. Let n2共r12兲 denote the number density of pairs of particles situated at distance r12 = 兩r1 − r2兩. Using the fact that n2共r12兲 = 0 for r12 ⬍ , we introduce a dimensionless function y 2共r12兲 defined by n2共r12兲 = n2共r12 − 兲y 2共r12兲,
共1兲
where is a unit step function. We assume that n2共r12兲 → n2 when r12 → ⬁. The twoparticle dimensionless correlation function h2共r12兲 is then defined by the cluster decomposition y 2共r12兲 = 1 + h2共r12兲, so that n2共r12兲 = 共r12 − 兲关1 + h2共r12兲兴, n2
共2兲
h2共r12兲 is thus supposed to satisfy the asymptotic condition lim h2共r12兲 = 0.
r12→⬁
共3兲
The second equilibrium YBG hierarchy equation establishes an exact relation between n2共r12兲 and the reduced three-particle number density n3共r1 , r2 , r3兲. Introducing the excluded volume factor we write the three-particle density as
133, 164507-1
© 2010 American Institute of Physics
164507-2
J. Chem. Phys. 133, 164507 共2010兲
Piasecki, Szymczak, and Kozak
n3共r1,r2,r3兲 = 共r12 − 兲共r13 − 兲 ⫻共r23 − 兲n3y 3共r12,r13,r23兲,
共4兲
The function n3 depends only on the distances rij = 兩rij兩 = 兩ri − r j兩, and is a symmetric function of the three variables r12 , r13 , r23. In the case of hard disks y 2 is related to y 3 through the second YBG hierarchy equation 共see Appendix A兲 d y 2共r兲 = n dr
冕
Our aim is to derive from Eq. 共10兲 the equivalent integral equation satisfied by the dimensionless correlation function
To this end, we insert Eq. 共11兲 into Eq. 共10兲 finding d ln关1 + H共x兲兴 dx
ˆ 共rˆ · ˆ 兲共兩r − ˆ 兩 − 兲 d
ˆ 兩兲, ⫻y 3共r, ,兩r −
= n2关1 + H共1兲兴
共5兲
ˆ and rˆ are unit vectors 兩 ˆ 兩 = 兩rˆ 兩 = 1. Putting rˆ · ˆ where = cos we rewrite Eq. 共5兲 in an explicit form d y 2共r兲 = n dr
冕
2
+
d cos y 3共r, , 冑r2 − 2r cos + 2兲 共6兲
When writing Eq. 共6兲, the equality ˆ 兩 − 兲 = 共r − 2 cos兲, 共兩r −
冕
冕
2
2
冎
⬁
dz共2 − z兲
x
冑 冉冊 冋 冑 冉冊 1−
= 共2 − x兲 −
x 2
where
x 2
2
z 2
2
1−
x 2
2
+ arccos
册
x , 2
共14兲
冋 冑 冉冊 x 2
1−
x 2
2
+ arccos
册
x , 2
共15兲
共8兲
Iⴱⴱ共x兲 = −
1 2
冕 冕 ⬁
2
dz
x
d
0
⫻cos 共z − 2 cos 兲H共冑z2 + 1 − 2z cos兲. 共16兲
共9兲
we find that it satisfies the nonlinear equation
冕
2
d cos Y共冑x2 − 2x cos + 1兲
It turns out that the angular integration in the formula for Iⴱⴱ共x兲 can be explicitly performed 共the calculation is presented in Appendix B兲. One eventually finds
0
⫻共x − 2 cos兲,
,
and
It is convenient to rewrite Eq. 共8兲 using the dimensionless distance x = r / . Denoting by Y共x兲, the function Y共x兲 = y 2共x兲,
1−
ln关1 + H共x兲兴 = 2n2关1 + H共1兲兴关Iⴱ共x兲 + Iⴱⴱ共x兲兴,
Iⴱ共x兲 = 共2 − x兲 −
⫻共r − 2cos兲.
冑 冉冊
we get
共7兲
dcos y 2共冑r2 − 2r cos + 2兲
共12兲
d共x − 2 cos兲cos = − 2共2 − x兲
0
d ln Y共x兲 = n2Y共1兲 dx
d cos H共冑x2 − 2x cos + 1兲
共13兲
d y 2共r兲 = ny 2共r兲y 2共兲 dr ⫻
2
0
III. THE KIRKWOOD SUPERPOSITION APPROXIMATION
Adopting Eq. 共7兲 leads to a closed equation
d cos 共x − 2 cos兲
0
In order to derive an integral equation for H共x兲 we integrate Eq. 共12兲 over the spatial interval 共x , ⬁兲. Using the formulae
冕
The Kirkwood superposition approximation consists in replacing in Eq. 共6兲 the three-particle density by the product of two-particle densities corresponding to three different pairs of particles
2
⫻共x − 2 cos兲 .
has been used. The rigorous relation 共6兲, which is valid for r ⬎ , will be the starting point for subsequent considerations.
ˆ 兩兲 → y 2共r兲y 2共兲y 2共兩r − ˆ 兩兲. y 3共r, ,兩r −
冕
再冕
0
0
⫻共r − 2 cos兲.
共11兲
H共x兲 = Y共x兲 − 1.
共10兲
valid in the region x ⱖ 1. Equation 共10兲 represents the closure of the YBG hierarchy corresponding to the superposition approximation.
Iⴱⴱ共x兲 = −
冕
x+1
x−1
ds共s − 1兲sH共s兲arccos
冉
冊
x2 + s2 − 1 . 2xs 共17兲
In this way, we arrive at an integral equation
164507-3
J. Chem. Phys. 133, 164507 共2010兲
Structural transition in the hard disk fluid
H
becomes slower as the surface fraction is increased, and a pronounced peak structure appears.
2
IV. LINEARIZATION OF KIRKWOOD’S EQUATION: THE RING APPROXIMATION 1
x 2
4
6
8
FIG. 1. Pair correlation function H共x兲 for the surface fraction = n2 / 4 = 0.35 共dashed line兲 and = 0.59 共solid line兲.
ln关1 + H共x兲兴 = 2n2关1 + H共1兲兴
再 冋
⫻ 共2 − x兲 − −
冕
x 2
冑 冉冊 x 2
1−
2
+ arccos
x 2
册
x+1
ds共s − 1兲sH共s兲
x−1
⫻arccos
冉
x2 + s2 − 1 2xs
冊冎
y 3共r12,r13,r23兲 = 1 + h2共r12兲 + h2共r13兲 + h2共r23兲
共19兲
where L is the integral operator given by L共H共x兲兲 ⬅ − 1 + exp 2n2关1 + H共1兲兴
再 冋
−
冕
冑 冉冊 x 1− 2
2
x + arccos 2
册
ds共s − 1兲sH共s兲
⫻arccos
冉
x2 + s2 − 1 2xs
冊冎 冊
y 3共r12,r13,r23兲 ⯝ 1 + h2共r12兲 + h2共r13兲 + h2共r23兲.
It is important here to note that while rejecting in the cluster decomposition 共22兲 the three-particle correlations we nevertheless retain in the hierarchy equation, the full excluded volume factor represented by the product of unit step functions 关see definition 共4兲兴. In fact, this factor represents the exact lowest order term in the density expansion of the threeparticle number density, and is of fundamental importance for the correct description of hard disks at low densities. Comparing Eq. 共23兲 with the superposition formula y 3共r12,r13,r23兲 ⯝ 关1 + h2共r12兲兴关1 + h2共r13兲兴关1 + h2共r23兲兴, 共24兲
.
共20兲
再
d H共x兲 = n2 − 2共2 − x兲 dx
The above integral equation for H共x兲 was solved by a standard Neumann method with successive over-relaxation.17 The iterative solutions are then given by Hn = 共1 − ␣兲Hn−1 + ␣L共Hn−1兲.
共23兲
we see that the ring approximation corresponds exactly to the linearization of the Kirkwood theory in h2. The linearized form of Eq. 共12兲 reads
x+1
x−1
共22兲
Neglecting h3 is thus equivalent to the approximation
H共x兲 = L共H共x兲兲,
x ⫻ 共2 − x兲 − 2
+ h3共r12,r13,r23兲.
共18兲
,
representing the superposition closure for the correlation function H共x兲. Our method of determining H共x兲 is based on the fact that the solution of Eq. 共18兲 satisfying the boundary condition limx→⬁ H共x兲 = 0 can be obtained numerically by iterations. Note that we can rewrite Eq. 共18兲 in the following form:
冉
Before continuing the analysis based on Eq. 共18兲, let us make a comment on the relationship between the superposition approximation and the ring approximation, well known from the kinetic theory 共see Ref. 15 and references given therein兲. Originally, the ring approximation was applied to the study of long wavelength hydrodynamic phenomena, and was defined by neglecting in the second equation of the dynamical BBGKY hierarchy all contributions from the threeparticle correlations. The three-particle correlation function h3 is defined by the cluster decomposition of the density y 3
+
2
1−
x 2
2
关1 + H共1兲 + H共x兲兴
dcos 共x − 2 cos兲
0
冎
⫻H共冑x2 + 1 − 2x cos兲 .
共21兲
The relaxation parameter ␣ was taken to be 0.25. The iterations were continued until successive values of H共x = 1兲 differed by less than ⑀ = 10−5, except in the vicinity of the threshold surface fraction ⴱ 共see Sec. V B兲, where the convergence was slow and the iterations were discontinued at ⑀ = 10−2. Examples of the correlation functions obtained in this way are presented in Fig. 1. As it is seen, the decay of H共x兲
冕
冑 冉冊
共25兲
In order to derive an integral equation for H共x兲, we proceed as before by integrating Eq. 共25兲 over the spatial interval 共x , ⬁兲. The result reads H共x兲 = 2n2关Jⴱ共x兲 + Jⴱⴱ共x兲兴, where
164507-4
J. Chem. Phys. 133, 164507 共2010兲
Piasecki, Szymczak, and Kozak
Z
Z共兲 =
10
p = 1 + n2共1 + H共1兲兲 = 1 + 2共1 + H共1兲兲, nkT 2 共30兲
where the contact value H共1兲 is also a function of . For comparison, we include here Z共兲 dependence as predicted by the scaled particle theory 共SPT兲18
5
0
0.1
ZSPT共兲 = x
0 0.2
0.3
0.4
0.5
0.6
FIG. 2. Compressibility factor Z vs the surface fraction for the full Kirkwood approximation 共circles兲, the ring approximation 共squares兲, and the scaled particle theory 关Eq. 共31兲, dashed line兴.
Jⴱ共x兲 = 共2 − x兲
再冋
冑 冉冊 册 冕 冑 冉冊 冎
x 2
−
2
1−
x 2
dz
1−
2
⫻关1 + H共1兲兴 +
x
+ arccos
x 2
关H共z兲兴 ,
共26兲
and Jⴱⴱ共x兲 = −
1 2
冕 冕 ⬁
2
dz
d cos 共z − 2 cos兲
0
x
⫻H共冑z2 + 1 − 2zcos兲 = Iⴱⴱ共x兲, 关see Eq. 共16兲兴. Using again the calculation presented in Appendix B, we eventually find Jⴱⴱ共x兲 = −
冕
x+1
ds共s − 1兲sH共s兲arccos
x−1
冉
冊
x2 + s2 − 1 . 2xs 共27兲
H共x兲 = 2n2共2 − x兲
再冋
冑 冉冊 冕 冑 冉冊 x 1− 2
2
⫻关1 + H共1兲兴 + − 2n2
冕
dz
x
1−
z 2
2
x + arccos 2
2
关H共z兲兴
x+1
x−1
ds共s − 1兲sH共s兲arccos
As it is seen, the iteration results agree with the SPT predictions up to approximately = 0.4.19 For larger packing fractions, the Kirkwood approximation tends to underestimate the compressibility factor, whereas the ring approximation overestimates it.
A. Breakdown of the method used for attractive interactions
The fundamental information concerning the internal structure of the system is contained in the spatial dependence of correlations. In particular, the law governing the asymptotic vanishing of correlations is of primary importance. In order to determine the behavior of H共x兲 at large distances, it seems natural to follow the moment analysis presented in Refs. 13 and 14. The calculation would proceed as follows. For x ⬎ 2, the Kirkwod equation 共12兲 can be conveniently written as d ln关1 + H共x兲兴 = n2关1 + H共1兲兴 dx
ˆ 共xˆ · ˆ 兲H共兩x − ˆ 兩兲. d
冉
2
冎
册
ˆ 兩兲 around When x Ⰷ 1, the power series expansion of H共兩x − the point 兩x兩 = x yields nonzero contributions only from terms ˆ 兲. The calculation up to involving odd powers of cos = 共xˆ · the third derivative of H yields the expansion ln关1 + H共x兲兴⬘ = n2关1 + H共1兲兴
2
冊
⫻
x +s −1 , 2xs
representing the ring approximation. Equation 共28兲 can again be solved numerically by iterations, and is expected to yield relevant results in the low density regime. The comparison of the results obtained with the full Kirkwood approximation and its linearization is presented in Fig. 2. Defining the surface fraction as
2 , 4
we plot here the compressibility factor Z共兲 given by
共29兲
冕
2
再
cos − cos H⬘共x兲
0
冉 冊
1 H⬘共x兲 ⬘ − 关cos − 共cos兲3兴 2 x
共28兲
= n
冕
共32兲
In this way we arrive at an integral equation x − 2
共31兲
V. ASYMPTOTIC DECAY OF CORRELATIONS: PREDICTING A STRUCTURAL TRANSITION
2
z 2
1 . 共1 − 兲2
冎
1 − 共cos兲3H共x兲 + ¯ , 6
共33兲
where ⬘ denotes the derivative with respect to x. Using the boundary condition limx→⬁ H共x兲 = 0 together with the equalities
冕
2
0
d共cos兲2 = ,
冕
2
0
3 d共cos兲4 = , 4
and adopting for large distances the asymptotic formula
164507-5
J. Chem. Phys. 133, 164507 共2010兲
Structural transition in the hard disk fluid
ln关1 + H共x兲兴 ⬃ H共x兲,
ReΚ,ImΚ 5
we find a linear differential equation of the form 1 H⬙共x兲 + H⬘共x兲 + ␣2H = 0, x
共34兲
A
0 4
where
冉
␣2 = 8 1 +
冊
1 ⬎ 0. n2关1 + H共1兲兴
B. Prediction of a structural transition
Let us consider again the region x ⬎ 2 where the Eq. 共32兲 holds. As limx→⬁ H共x兲 = 0, we can replace in Eq. 共32兲 the function ln关1 + H共x兲兴 by H共x兲, and consider the equation
冕
2
d cos H共冑x2 − 2x cos + 1兲,
12
16
20
5
However, Eq. 共34兲 is the equation for the Bessel function J0共␣x兲. The solution of Eq. 共34兲 is thus an oscillating function whose amplitude decays as 1 / 冑x. Unfortunately, the above result is in disagreement with numerical predictions of exponentially damped oscillations. Use of the moment expansion as developed in Ref. 13 leads to erroneous results. Clearly one has to take into account the ˆ 兩兲 to get whole infinite series in the expansion of H共兩x − reliable predictions. We have thus to give up this kind of expansion, and look for a different approach. In a broad outline, the failure of the moment expansion developed in Refs. 13 and 14 to describe the structural transition in the hard disk fluid can be traced to the nature of the governing intermolecular potential. The studies13,14 deal with the description of correlations in the vicinity of the liquidvapor critical point, where attractive forces play a dominant role. Using the moment expansion, one recovers at large distances the classical Ornstein–Zernike formula. The present study deals with the fusion transition, where short-range, repulsive forces play the critical role. The new method introduced here 共see Sec. V B兲 effectively accounts for the difference in the potential governing these two transitions and, anticipating our later results, leads to an analytic prediction of the packing fraction at which a structural transition occurs in the hard disk 共and hard sphere兲 fluid.
A d H共x兲 = dx 2
8
共35兲
0
FIG. 3. The real part 共solid line兲 and imaginary part 共dashed line兲 of the root of Eq. 共38兲 corresponding to the slowest decaying mode.
notice that all derivatives of H共x兲 satisfy the same equation. It is thus natural to consider H共x兲 as a linear combination of exponential modes exp共x兲, where is a complex number. The function exp共x兲 satisfies Eq. 共37兲 provided solves the equation
=
A 2
冕
2
d cos exp共− cos兲 = − AI1共兲,
共38兲
0
where I1共兲 is a modified Bessel function. The physically acceptable solutions = a + ib are those with a negative real part a ⬍ 0 that assures the exponential damping of oscillations. The first root 共with the smallest absolute value of a兲 corresponding to the slowest decay of the correlation function is shown in Fig. 3. The values of 共兲 predicted with the use of Eq. 共38兲 are in good agreement with the decay of the amplitude of H共x兲 determined by the iterative solution of integral equation 共19兲. Plotting the absolute value of H共x兲 on a logarithmic plot, and fitting it by the single mode ␣e⬘x, we obtain values of ⬘ that are slightly below those obtained by solving Eq. 共38兲. For example, for = 0.3 we get ⬘ ⬇ −2.1 共cf. Fig. 4兲, whereas the corresponding value of for that surface fraction is ⬇ −2.02. Similarly, for = 0.55, we get, respectively, ⬘ ⬇ −0.8 and ⬇ −0.7. The fact that the values of ⬘ remain slightly below those of can be understood by noting that corresponds to the slowest decaying mode, whereas in the numerical data on H共x兲 we also see nonzero contributions from other, faster decaying modes. An interesting feature of the 共兲 dependence presented in Fig. 3 is the fact that the root becomes purely imaginary at
where H 1
A = 2n2关1 + H共1兲兴. We then use the expansion 2 冑x2 − 2x cos + 1 = x − cos + sin + ¯ , 2x
0.01
共36兲
106
to arrive at the equation A d H共x兲 = dx 2
冕
2
d cos H共x − cos兲,
104
108
共37兲
1010
0
valid for x Ⰷ 1. In order to determine the large x behavior of correlations, we have thus to analyze the solution of Eq. 共37兲. We
2
4
6
8
10
x
FIG. 4. The absolute value of H共x兲 for = 0.3 together with a corresponding fit of the form ae⬘x with ⬘ = −2.1.
164507-6
J. Chem. Phys. 133, 164507 共2010兲
Piasecki, Szymczak, and Kozak
8 ⬎ 15.12,
A 16
13
12
共41兲
representing for a physically impossible condition 共beyond close packing兲. The ring approximation is unable to describe a qualitative change in the hard disk correlation function. For any accessible surface fraction, it predicts exponentially damped oscillations.
15 14
⬎ 1.89,
or
11
10
C. Structural transition in a hard sphere fluid
x 0.5
0.52
0.54
0.56
0.58
0.6
0.62
FIG. 5. The extrapolation of A共兲 dependence. The points correspond to the values of A obtained from the iterative procedure, and the dashed line is given by A = Aⴱ ⬇ 15.12.
A = Aⴱ ⬇ 15.1. This point can be made more precise by analyzing when Eq. 共38兲 acquires a purely imaginary solution = ib. As I1共ib兲 = iJ1共b兲, Eq. 共38兲 implies the condition 1 J1共b兲 1 = 关J0共b兲 + J2共b兲兴 = − , 2 b A
共39兲
where the first equality follows from the recurrence relation between Bessel functions Jn. The absolute minimum of the sum of Bessel functions 关J0共b兲 + J2共b兲兴 equals ⫺0.1323. The necessary condition for the disappearance of damping has thus the form 2 ⱕ 0.1323, A
or
2 = 4n3
冉
共40兲
The above inequality shows that the nature of correlations could change for sufficiently high values of the surface fraction 关Eq. 共29兲兴 occupied by hard disks. The determination of the value of ⴱ corresponding to the equality Aⴱ = 8ⴱ关1 + H共x = 1; ⴱ兲兴 = 15.12, requires the knowledge of the contact value H共1兲 as a function of the surface fraction. We have studied this question numerically. An accurate estimation of the precise value of ⴱ corresponding to Aⴱ is difficult because as we approach ⴱ, the iteration procedure demands an increasingly larger number of iterations to converge to a solution with the required accuracy. Additionally, a computational domain over which the solution is sought must also be progressively extended as we approach ⴱ, since the decay of H共x兲 is very weak there. We estimated the value of ⴱ by calculating A共兲 for several values of in the range of 0.5ⱕ ⱕ 0.6 and then extrapolating to larger values of , as illustrated in Fig. 5. In this way, we obtain the estimate of ⴱ ⬇ 0.622. Let us close this section with a comment on the ring approximation. In order to find the asymptotic behavior of correlations satisfying Eq. 共28兲, it is sufficient to consider the linearized version of Eq. 共35兲. However, this is equivalent to the replacement of the factor A = 2n2关1 + H共1兲兴 by 2n2 = 8. The inequality 共40兲 becomes then
冊
sh共兲 − ch共兲 .
共42兲
In passing to the Kirkwood superposition approximation, we need only to replace the factor 4n3 in the above equation by = 4n3关1 + H共1兲兴. This fact follows from the previously made remark that the ring approximation represents exactly the linearized Kirkwood theory. In order to determine the range of values of ⬎ 0 beyond which the exponential damping becomes impossible we explore the possibility of vanishing of the real part a of the complex number = a + ib. Equation 共42兲 reduces then to
冉
b2 = cosb −
A = 2n2关1 + H共1兲兴 ⱖ Aⴱ
Aⴱ = 15.1171 ¯ ⯝ 15.12.
The integral equation satisfied by the correlation function H共x兲 of a three-dimensional hard sphere fluid for x Ⰷ 1 has been analyzed in Ref. 15 within the ring approximation. The frequencies of the exponential modes exp共x兲 describing the long distance decay of correlations are in this case solutions of the equation 关see Eq. 共22兲 in Ref. 15兴
冊
sin b , b
or
=
b3 . b cosb − sin b
共43兲
The absolute minimum of the function y共b兲 =
b3 , b cosb − sin b
共44兲
in the region where y共b兲 ⬎ 0 equals 共45兲
y min = 34.81 ¯ .
Hence, for ⬎ 34.81, Eq. 共42兲 acquires purely imaginary solutions and the exponential damping vanishes. It is quite remarkable that 60 years ago, Kirkwood et al.6 concluded from numerical studies of the integral equation for H共x兲 that when exceeds 34.8 the correlation function H共x兲 is not integrable anymore. Our method provides a simple analytic confirmation of this result.
D. The hard rod fluid: Testing the method
The rigorous calculation of the two-particle correlation function for hard rods 共see e.g., Ref. 20兲 shows that its structure corresponds to exponentially damped oscillations at all possible densities. One finds ⬁
n关H共x兲 + 1兴 = 兺 关x − 共k + 1兲兴 k=0
k关x − 共k + 1兲兴k k!
⫻exp兵− 关x − 共k + 1兲兴其, where
共46兲
164507-7
J. Chem. Phys. 133, 164507 共2010兲
Structural transition in the hard disk fluid
Κ ns
0 0.2
0.4
0.6
0.8
1
H 5
2
4
4
3 2
6
1
x
8
2
FIG. 6. The real part of n = − + Wn共 exp 兲, n = 1 , . . . , 5 共top to bottom兲 as a function of n = / 共1 + 兲.
=
n . 1 − n
In fact, the superposition law turns out to be exact for a one-dimensional hard rod fluid,20 and the second equation of the equilibrium hierarchy takes a particularly simple form H⬘共x兲 = n关共x − 2兲H共x − 1兲 − H共x兲 − 共2 − x兲兴 ⫻关H共1兲 + 1兴.
共47兲
The formula 共46兲 represents the solution of Eq. 共47兲. When x ⬎ 2, we find H⬘共x兲 = n关H共x − 1兲 − H共x兲兴关H共1兲 + 1兴.
共49兲
It turns out that all solutions of Eq. 共49兲 can be expressed in terms of the multivalued Lambert W function. Indeed, the special function W共z兲 is defined on the complex plane by the equation z = W共z兲exp关W共z兲兴.
共50兲
But, Eq. 共49兲 can be rewritten as 共 + 兲exp共 + 兲 = exp .
共51兲
It follows that
= − + W共 exp 兲.
6
8
0 FIG. 7. The exact form of the correlation function for hard-rod fluid 共solid兲 at n = 0.8 and its asymptotic form H共x兲 = Aex with 共兲 calculated using Eq. 共52兲.
As illustrated in Figs. 7 and 8, the correlation function H共x兲 rapidly approaches the asymptotic form given by the slowest decaying mode Ae1x. Note that not only the exponential decay rate but also the oscillation period of the function agree with that given by Ae1x starting from x ⬇ 3. This shows that the other modes play a negligible role in influencing the behavior of H共x兲 for intermediate and large x values, thus lending further support to our approach of focusing on the slowest decay mode only.
共48兲
Applying the same method as that used for hard disks, we look for exponential modes exp共x兲 solving Eq. 共48兲. The complex frequency satisfies the equation
= n关H共1兲 + 1兴共e− − 1兲 = 共e− − 1兲.
4
共52兲
The principal branch of Lambert function W0共z兲 obeys W0共xex兲 = x for real x, thus = 0 for that branch. However, other branches Wn共z兲 with n ⬎ 0 give values of n with a negative real part corresponding to the decay of the correlation function. The first five solutions, n共兲 , n = 1 , . . . , 5, are presented in Fig. 6 as functions of n = / 共1 + 兲 共modes with larger n decay faster兲. The negative real part of n , n ⬎ 0 is different from zero for any accessible density, and vanishes only at close packing where n = 1. There are thus no purely imaginary solutions = ib of Eq. 共49兲. This confirms the correctness of the method, showing the absence of any structural transition in one dimension over the whole range of densities below close packing.
VI. DISCUSSION AND CONCLUSIONS
The question of whether a system of particles interacting via a purely repulsive potential 共only兲 can undergo a phase transition has been under continuous investigation since first posed and addressed by Kirkwood over 70 years ago. For a system of hard disks, the first numerical evidence was provided by Alder and co-workers.21,22 The data reported in Ref. 21 showed that the hard disk freezing transition occurred at a density smaller than the density of closest packing 关corresponding to an area fraction of 0 = / 冑12= 0.906 90兴, and suggested that the liquid to solid transition was first order. The most recent Monte Carlo simulations 共on a system of 4 ⫻ 106 disks兲 of Mak16 suggest that melting consists of a continuous transition from the ordered solid to an intermediate 共hexatic兲 phase23–27 at a packing fraction = 0.723, and
H 1 0.1 0.01 0.001 104 105 2
4
6
8
10
12
14
x
FIG. 8. The absolute value of H共x兲 − 1 for hard rod fluid 共solid兲 at n = 0.8 and the exponential asymptote H共x兲 = AeRe共兲x.
164507-8
J. Chem. Phys. 133, 164507 共2010兲
Piasecki, Szymczak, and Kozak
either a very weak first-order or a continuous transition from the intermediate phase to the fluid phase at a packing fraction = 0.699. From the analysis and numerical evidence presented in Sec. V B, we have calculated the area fraction at which a structural change in the hard disk fluid can take place, viz., ⴱ ⬃ 0.622. Converting this area fraction to a packing fraction gives ⴱ = ⴱ / 0 ⬃ 0.718. When compared to the estimate reported by Mak for the transition from the ordered solid phase to a dense fluid phase, ⬃ 0.723, one finds the two values are in substantial agreement. Also of interest is our prediction of a structural transition in a system of hard spheres. As noted in Sec. V C, Kirkwood, Maum and Alder6 found that for values of ⱖ 34.8, no solutions of the YBG and Kirkwood integral equations exist for which, in their notation, x2关g共x兲 − 1兴 is integrable. We find that the value of beyond which exponential damping becomes impossible is = 34.81. Hence, the results for hard spheres appear to be in exact agreement. When we consider a system of hard rods, the analytic method developed here confirms the absence of any structural transition in d = 1 for the entire range of densities below close packing 关a result that was already known to Rayleigh 共Ref. 28兲 and Korteweg 共Ref. 29兲兴. Our prediction of a structural phase transition is based on the analysis of an integral equation whose derivation assumed sufficiently fast decay of correlations. If a phase becomes ordered and correlations do not decay, the integral equation to which our method was applied does not hold. This can preclude the possibility of studying a region where a new equilibrium phase may be formed. In particular, the identification of a region intermediate between melting and freezing 共see following paragraph兲 and the characterization of an ordered 共solid兲 phase is certainly beyond the scope of our approach. To study these questions within the Kirkwood superposition approximation, one must go back to the original YBG hierarchy equation 共Kirkwood’s closure does not assume the rapid decay of correlations兲. To develop this point further, there are three structural aspects of the hard disk transition that are not captured by the method developed in this paper. First, we find no evidence for the existence of an intermediate, or hexatic, phase predicted by the Kosterlitz, Thouless, Halperin, Nelson, and Young 共KTHNY兲 theory of d = 2 melting,23–27 and supported by Mak’s Monte Carlo simulations. Second, we find no evidence for the development of a shoulder on the second maximum of disk radial distribution function in the vicinity of the freezing transition 共 = 0.686兲 reported by Truskett et al.30 based on molecular dynamics simulations, and later correlated with structural rearrangements occurring at increasing disk density.31–34 Third, we cannot confirm the existence of regions of fivefold coordination in the dense fluid phase, first predicted by Bernal35–37 based on his “ball and spoke” model of a random assembly of hard-core particles, although it has been conjectured that the hexatic phase might be correlated with randomly dispersed regions of 5-, 6-, and 7-member disk clusters forming percolated tessellations that span the d = 2 space.38 The larger point, however, relates to the original Kirkwood prediction, viz., that a system of particles interacting
via purely repulsive forces, here hard disks but also hard spheres, can undergo a phase transition. Although not widely accepted at first, following the work of Onsager on the isotropic-nematic transition in a d = 3 dimensional system of thin hard rods,39 there developed a gradual realization that a phase transition can be entropy driven. As elaborated by Frenkel,40 in hard-core systems the entropy in the ordered crystalline phase is actually larger than in the fluid phase; quoting directly, “the entropy decreases because the density is no longer uniform in orientation or position, but the entropy increases because the free-volume per particle is larger in the ordered than in the disordered phase.” The present contribution provides further evidence for the essential correctness of Kirkwood’s insight.
APPENDIX A: DERIVATION OF EQ. „5… FROM THE BBGKY HIERARCHY
Consider a gas of hard disks of mass m and diameter . We denote by j ⬅ 共r j , v j兲, j = 1 , 2 , . . ., the one-particle state in which a disk has position r j and velocity v j. The average number density of s-particle clusters occupying at time t the s-particle state 共1 , 2 , . . . , s兲 is called the s-particle reduced distribution f s共1 , 2 , . . . , s ; t兲. The dynamical evolution of the hard disk fluid is described in the thermodynamic limit by the BBGKY hierarchy equations. The second of them establishes a relation between f 2 and f 3
冉
冊
+ v1 · + v2 · − ¯T共1,2兲 f 2共1,2;t兲 t r1 r2 =
冕
¯ 共1,3兲 + ¯T共2,3兲兲f 共1,2,3;t兲. d3共T 3
共A1兲
The effects of binary collisions are described by the operator ¯T共i , j兲 ¯T共i, j兲 =
冕
ˆ v12 · ˆ 共vij · ˆ 兲关␦共rij − 兲bˆ d
− ␦共rij + 兲兴.
共A2兲
Here the Dirac ␦-distributions restrict the distances 兩rij兩 between the centers of the disks at the moment of impact to ˆ is oriented perpentheir diameter = 兩兩. The vector = dicularly to the surface of colliding disks at the point of impact. The action of the operator bˆ consists in replacing the velocities vi , v j by their precollisional values v⬘i , v⬘j corresponding to the inverse elastic collision. As in elastic collisions the kinetic energy is conserved, in the case of products of Maxwell distributions
共v兲 =
冉
冊
m exp共− mv2/2kBT兲, 2 k BT
we find bˆ 关共vi兲共v j兲兴 = 共v⬘i 兲共v⬘j 兲 = 共vi兲共v j兲. Hence, in the case of equilibrium reduced distributions
共A3兲
164507-9
J. Chem. Phys. 133, 164507 共2010兲
Structural transition in the hard disk fluid
f s共1,2, . . . ,s兲 = ns共r1,r2 ¯ rs兲共v1兲共v2兲 ¯ 共vs兲, 共A4兲 Eq. 共A1兲 takes the form
冉
v1 · + v2 · − r1 r2
冕
ˆ v12 · ˆ ␦共rij − 兲 d
冊
冕
dr3dv3
冕
− r12
冕
冕 冕 dr3
共A5兲
冊
Iⴱⴱ共x兲 = −
冕
ˆ v12 · ˆ ␦共r12 − 兲, d
y 2共r12兲. r12
冕
x−cos
ds共s − cos兲 共B2兲
d
ds共x + 冑1 − 2 − s兲共s − x + 冑1 − 2兲
⫻H共冑s2 + 2兲共s − 冑1 − 2兲.
共A7兲
共B3兲
We now introduce a new integration variable
共A6兲
z = 冑s2 + 2 . As zdz = sds, we find
冕 冕 1
H共z兲共冑z2 − 2 2 冑 z − 0 − 冑1 − 2兲共冑1 − 2 − 兩x − 冑z2 − 2兩兲
Iⴱⴱ共x兲 = − 共A8兲
=−
d
冕 冕
dz
z
2
dz zH共z兲共z − 1兲 1
d
⫻
1
冑z2 − 2 共
冑1 − 2 − 兩x − 冑z2 − 2兩兲.
ˆ v12 · ˆ 共兩共r12 − 兲兩 − 兲 d
共B4兲 共A9兲
As the equality L = R must hold for any value of the relative velocity v12, we finally find 共when r12 ⱖ 兲
冕
冕 冕 1
0
⫻y 3共r12, ,兩共r12 − 兲兩兲.
dy 2共r12兲 = n dr12
x+cos
0
On the right-hand side owing to the presence of ␦-distributions, we can perform integration over variable r3, thus obtaining R = n共r12 − 兲
冕
d cos
Putting then = sin , we get
which reduces the left hand side of Eq. 共A6兲 to L = 共r12 − 兲v12 ·
/2
⫻H共冑s2 + sin2 兲.
We now use the identity
共r12 − 兲 ⬅ r12
冕
0
⫻共r23 − 兲 − 兲y 3共r12, ,r23兲
v12 ·
dcos 关共z − 2 cos兲
0
Changing the order of integrations with the use of the asymptotic decay of the correlation function we arrive at a convenient formula
Iⴱⴱ共x兲 = −
ˆ 关共v1 · ˆ ␦共r13 − 兲 d
ˆ ␦共r23 − 兲共r13 − 兲y 3共r12,r13, 兲兴. + v2 ·
/2
dz
共B1兲
ˆ v12 · ˆ ␦共r12 − 兲 共r12 − 兲y 2共r12兲 d
= n共r12 − 兲
冕 冕
⫻H共冑z2 + 1 − 2z cos兲 − H共冑z2 + 1 + 2z cos兲兴.
ˆ 关v13 · ˆ ␦共r13 − 兲 d
We can divide both sides of Eq. 共A5兲 by the product 共v1兲共v2兲, and perform the integration over the v3 variable. Moreover, we introduce explicitly the excluded volume factors by using Eqs. 共1兲 and 共4兲. The hierarchy equation becomes v12 ·
d cos 共z − 2 cos兲
0
x
x
⫻共v1兲共v2兲共v3兲.
2
dz
⬁
=−
ˆ ␦共r23 − 兲兴n3共r12,r13,r23兲 + v23 ·
冉
冕 冕 ⬁
⫻H共冑z2 + 1 − 2z cos兲
⫻n2共r12兲共v1兲共v2兲 =
1 2
Iⴱⴱ共x兲 = −
ˆ rˆ 12 · ˆ 共兩r12 − 兩 − 兲 d
⫻y 3共r12, ,兩r12 − 兩兲,
Here the inequality 冑1 − 2 ⬎ 兩x − 冑z2 − 2兩 is equivalent to
冑z2 − 2 ⬎ x
2
+ z2 − 1 , 2x
which leads to the formula 共A10兲
which is Eq. 共5兲 of Sec. II.
Iⴱⴱ共x兲 = −
冕 冕
dzzH共z兲共z − 1兲 1
⫻ APPENDIX B: DERIVATION OF EQ. „17…
We perform here the angular integration in the contribution to the correlation function
0
d
1
冑z2 −
2
冉冑
z2 − 2 −
冊
x2 + z2 − 1 . 2x 共B5兲
Putting = z, we find
164507-10
J. Chem. Phys. 133, 164507 共2010兲
Piasecki, Szymczak, and Kozak
Iⴱⴱ共x兲 = −
冕 冕
1/z
d
⫻
0
1
冑1 −
2
9
冉冑
1 − 2 −
冊
x2 + z2 − 1 . 共B6兲 2xz
The change of the integration variable w = 冑1 − 2 yields the formula Iⴱⴱ共x兲 = −
冕 冕 冕 冕 冕
⫻ =−
=−
⫻
dz zH共z兲共z − 1兲 1
dw
冑1−1/z2 冑1 − w
⫻
冉
w− 2
x2 + z2 − 1 2xz
1
dw
冑1 − w
冉
1− 2
x2 + z2 − 1 2xz
冉
x2 + z2 − 1 − arcsin 2 2xz
冕
冊
dz zH共z兲共z − 1兲关1 − 共x − z兲2兴
冊冎
x+1
x−1
共B7兲
.
Finally, we arrive at the result Iⴱⴱ共x兲 = −
冊
dz zH共z兲共z − 1兲
共x2+z2−1兲/2xz
再
dz zY共z兲共z − 1兲arccos
冉
冊
x2 + z2 − 1 , 2xz 共B8兲
used in Eq. 共17兲. T. L. Hill, Statistical Mechanics 共McGraw-Hill, New York, 1956兲. S.A. Rice and P. Gray, The Statistical Mechanics of Simple Liquids 共Interscience, New York, 1965兲. 3 J. G. Kirkwood, J. Chem. Phys. 3, 300 共1935兲. 4 J. G. Kirkwood and E. Monroe, J. Chem. Phys. 9, 514 共1941兲. 5 J. G. Kirkwood and E. M. Boggs, J. Chem. Phys. 10, 394 共1942兲. 6 J. G. Kirkwood, E. K. Maum, and B. J. Alder, J. Chem. Phys. 18, 1040 共1950兲. 7 J. P. Hansen and I. R. McDonald, Theory of Simple Liquids, 3rd ed. 共Academic, London, 2006兲. 1 2
J. D. Weeks, S. A. Rice, and J. J. Kozak, J. Chem. Phys. 52, 2416 共1970兲. J. J. Kozak, Adv. Chem. Phys. 40, 229 共1979兲. 10 B. Bagchi, C. Cerjan, and S. A. Rice, Phys. Rev. B 28, 6411 共1983兲. 11 A. D. J. Haymet, Annu. Rev. Phys. Chem. 38, 89 共1987兲. 12 M. Baus and C. F. Tejero, Equilibrium Statistical Mechanics: Phases of Matter and Phase Transitions 共Springer, New York, 2008兲. 13 G. L. Jones, J. J. Kozak, E. Lee, S. Fishman, and M. E. Fisher, Phys. Rev. Lett. 46, 795 共1981兲; M. E. Fisher and S. Fishman, ibid. 47, 421 共1981兲. 14 M. E. Fisher and S. Fishman, J. Chem. Phys. 78, 4227 共1983兲 共This reference provides an extensive account of the primarily analytic studies of the YBG equation using the moment expansion.兲; G. L. Jones, E. K. Lee, and J. J. Kozak, ibid. 79, 459 共1983兲 共This reference provides a similar account of the primarily numerical studies of the YBG equation.兲. 15 J. Piasecki and R. Soto, Physica A 379, 409 共2007兲. 16 C. H. Mak, Phys. Rev. E 73, 065104 共2006兲. 17 R. R. Kleinman and P. M. van den Berg, PIER 5, 67 共1991兲. 18 H. Reiss, H. L. Frisch, and J. L. Lebowitz, J. Chem. Phys. 31, 369 共1959兲. 19 K. Helfand, H. L. Frisch, and J. L. Lebowitz, J. Chem. Phys. 34, 1037 共1961兲. 20 Z. Salsburg, R. Zwanzig, and J. G. Kirkwood, J. Chem. Phys. 21, 1098 共1953兲. 21 B. J. Alder and T. E. Wainwright, Phys. Rev. 127, 359 共1962兲. 22 B. J. Alder, W. G. Hoover, and T. E. Wainwright, Phys. Rev. Lett. 11, 241 共1963兲. 23 J. M. Kosterlitz and D. J. Thouless, J. Phys. 共Paris兲 C5, 1124 共1972兲. 24 B. I. Halperin and D. R. Nelson, Phys. Rev. Lett. 41, 121 共1978兲. 25 D. R. Nelson and B. I. Halperin, Phys. Rev. B 19, 2457 共1979兲. 26 A. P. Young, Phys. Rev. B 19, 1855 共1979兲. 27 For a comprehensive review, see: D. R. Nelson, in Phase Transitions and Critical Phenomena, edited by C. Domb and J. L. Lebowitz 共Academic, London, 1983兲, Vol. 7. 28 Lord Rayleigh, Nature 共London兲 45, 80 共1891兲. 29 D. T. Korteweg, Nature 共London兲 45, 152 共1891兲. 30 T. M. Truskett, S. Torquato, S. Sastry, P. G. Debenedetti, and F. H. Stillinger, Phys. Rev. E 58, 3083 共1998兲. 31 A. Huerta, G. G. Naumis, D. T. Wasan, D. Henderson, and A. Trokhymchuk, J. Chem. Phys. 120, 1506 共2004兲. 32 F. Moucka and I. Nezbeda, Phys. Rev. Lett. 94, 040601 共2005兲. 33 A. Huerta and G. G. Naumis, Phys. Rev. Lett. 90, 145701 共2003兲. 34 A. Huerta, D. Henderson, and A. Trokhymchuk, Phys. Rev. E 74, 061106 共2006兲. 35 J. D. Bernal, Trans. Faraday Soc. 33, 27 共1937兲. 36 J. D. Bernal, Nature 共London兲 183, 141 共1959兲. 37 J. D. Bernal and S. V. King, Physics of Simple Fluids 共John Wiley & Sons, Inc., New York, 1968兲, p. 231. 38 J. J. Kozak, J. Brzezinski, and S. A. Rice, J. Phys. Chem. B 112, 16059 共2008兲. 39 L. Onsager, Ann. N.Y. Acad. Sci. 51, 627 共1949兲. 40 D. Frenkel, Physica A 263, 26 共1999兲. 8
dz zH共z兲共z − 1兲
