Self-Diffusion Coefficient of the Hard-Sphere Fluid: System Size

of the nonlinear least-squares Marquardt−Levenberg algorithm in the gnuplot ...... (44). Rosenfeld, Y. J. Phys. Condens. Matt. 1999, 11, 5415. [...
3 downloads 0 Views 203KB Size
Download PDF
J. Phys. Chem. B 2007, 111, 1455-1464

1455

Self-Diffusion Coefficient of the Hard-Sphere Fluid: System Size Dependence and Empirical Correlations D. M. Heyes* and M. J. Cass DiVision of Chemistry, School of Biomedical and Molecular Sciences, UniVersity of Surrey, Guildford GU2 7XH, U.K.

J. G. Powles School of Engineering, UniVersity of Wales Swansea, Singleton Park, Swansea SA2 8PP, U.K.

W. A. B. Evans The Physics Laboratory, UniVersity of Kent, Canterbury, Kent CT2 7NZ, U.K. ReceiVed: NoVember 8, 2006; In Final Form: December 18, 2006

Molecular dynamics simulations have been used to calculate the self-diffusion coefficient, D, of the hard sphere fluid over a wide density range and for different numbers of particles, N, between 32 and 10 976. These data are fitted to the relationship D ) D∞ - AN-R where the parameters D∞, A, and R are all densitydependent (the temperature dependence of D can be trivially scaled out in all cases). The value R ) 1/3 has been predicted on the basis of hydrodynamic arguments. In the studied system size range, the best value of R is ∼1/3 at intermediate packing fractions of ∼0.35, but increases in the low- and high-density extremes. At high density, the scaling follows more closely that of the thermodynamic properties, that is, with an exponent of order unity. At low packing fractions (less than ∼0.1), the exponent increases again, appearing to approach a limiting value of unity in the zero-density limit. The origin of this strong N dependence at low density probably lies in the divergence in the mean path between collisions, as compared with the dimensions of the simulation cell. A new simple analytical fit formula based on fitting to previous simulation data is proposed for the density dependence of the shear viscosity. The Stokes-Einstein relationship and the dependence of D on the excess entropy were also explored. The product Dηsp with p ) 0.975 was found to be approximately constant, with a value of 0.15 in the packing fraction range between 0.2 and 0.5.

1. Introduction Hard spheres have proved a useful model fluid for calculating the thermophysical and dynamical properties of real molecular fluids because their properties are largely determined by the repulsive part of the intermolecular potential for which the hard sphere interaction is the simplest representation (see, e.g., ref 1, Chapter 6). The hard sphere potential, φ(r), is

φ(r) )

{

∞ reσ 0 r>σ

(1)

where r is the interparticle separation and σ is the hard sphere diameter. The packing fraction, ζ, is a convenient measure of the density and is defined as ζ ) πNσ3/6V, where N is the number of hard spheres in volume V. The self-diffusion coefficient, D, is key to understanding the dynamics and relaxation processes in liquids, and hence, so is the self-diffusion coefficient of the HS fluid. Enskog derived a kinetic theory for D applicable at finite densities (see, e.g., ref 2). The first molecular dynamics, MD, calculations of the self-diffusion coefficient of the hard sphere fluid were reported in the pioneering paper by Alder and Wainwright.3 Corrections to the Enskog result derived from the MD simulations were proposed * To whom correspondence should be addressed. E-mail: d.heyes@ surrey.ac.uk.

by Alder and co-workers.5,4 There have been a number of studies since; for example, Erpenbeck and Wood,6 Woodcock,7 and a useful publication by Speedy, who collected the previous simulation data together.8 The latter two papers include data for D in the metastable fluid region, that is, at a packing fraction >0.493.9 A review of various simulation studies of the transport coefficients of hard spheres and their relationship to theory was given by Sigurgeirsson and Heyes.10 For a discussion of the relationship between the self-diffusion coefficients of hard spheres and other model fluids, see, for example, ref 11. To make meaningful comparisons with theory, we require the value of the self-diffusion coefficient in the thermodynamic limit. It has been known for some time that there is a significant system size dependence of the self-diffusion coefficient of the hard sphere system (as compared with the other transport coefficients). Finite size effects for a kinetic lattice-gas model were carried out by Kob and Andersen.12 Many of the earlier hard sphere MD studies were carried out on relatively small system sizes, typically on the order of N ) 500 particles (apart from those of Erpenbeck and Wood, who used up to N ) 4000 hard spheres.6) It is not practical to carry out simulations with a large enough N to essentially achieve the thermodynamic limit, and it is necessary to extrapolate D obtained for selected system sizes to the thermodynamic limit (i.e., N f ∞) for use in statistical mechanical models. Erpenbeck and Wood carried out an extensive theoretical and molecular dynamics study of the

10.1021/jp067373s CCC: $37.00 © 2007 American Chemical Society Published on Web 01/24/2007

1456 J. Phys. Chem. B, Vol. 111, No. 6, 2007

Heyes et al.

hard sphere velocity autocorrelation function, VACF, which is used in a Green-Kubo formula to obtain D (e.g., see ref 1, p 201), and the self-diffusion itself, from the dilute gas to near solidification.6 They considered the various contributions to the VACF and used this to extrapolate their finite system size simulations of the self-diffusion coefficient to the thermodynamic limit. The actual form of the N dependence was not explicitly considered in the Erpenbeck and Wood study, however, although it will, of course, be implicit in their results. Dunweg and Kremer13 and, more recently, Fushiki14 and Yeh and Hummer15 derived an expression for the N dependence of the self-diffusion coefficient, taking into account the effects of the periodic images using a hydrodynamics model, which gives an N-1/3 dependence for D(N). They also carried out MD simulations of hard spheres14 and Lennard-Jones and model water molecules15 at liquidlike densities. The present study follows from this work by focusing on hard sphere self-diffusion over a wide range of packing fractions. D(N, ζ) has been obtained for most of the fluid range, including very low densities, and these data are parametrized with relatively simple semiempirical analytic formulas. Apart from providing another estimate of D in the thermodynamic limit, the N dependence is examined in more detail, admittedly at a largely empirical level. The simulations were carried out for many system sizes, from N ) 32 to 10 976 particles, for typically 50 million or more hard sphere collisions. The self-diffusion coefficients were obtained from both the velocity autocorrelation function (VACF) and the mean square displacement (MSD) routes (e.g., see ref 16, p 95, and ref 17). Erpenbeck and Wood6 used a center of mass frame-corrected definition of the velocity for the VACF. In our study, the VACFs were determined in the laboratory frame, and center of mass corrections were accounted for in the formula used to describe the system size dependence of the self-diffusion coefficient. We confirmed that the laboratory frame VACFs decayed to 0 before truncation within the small statistical fluctuations in the VACF. We found the VACF route to give better statistics than the MSD method. The VACF formula uses the particle velocity Vi for particle of index i,

Figure 1. Plot of D(ζ, N) against N-1/3 for a set of MD diffusion coefficient values. The packing fraction, ζ ) πNσ3/6V, where σ is the hard sphere diameter and N is the number of hard spheres in volume V, from top to bottom is 0.25, 0.30, 0.35, 0.40, 0.45, 0.49, 0.50, and 0.51. The lines are linear regression least-squares fits to the data using eq 6 with R ) 1/3. The data span from N ) 10 976 with an abscissa value of 0.045 to N ) 32 with an abscissa value of 0.315.

(2)

Figure 2. As for Figure 1, except that the diffusion data is plotted against N-1/2 and R ) 1/2 in eq 6.

where 〈‚‚‚〉s denotes a time correlation or VACF. For practical reasons, the VACF cutoff time, tc, in eq 2 is finite due to computer resource limitations. As confirmed in ref 6, it is important to include the effects of any long-time tail in the VACF. Estimated contributions from hydrodynamic tails were not included in the diffusion coefficients because a sufficiently large value of tc was used. There was an insignificant contribution to the computed self-diffusion coefficient from a long-time tail for t > tc. The value of tc was taken to be at least several hundred collisions per particle, and for the highest densities, tc = 100 hard sphere time units. The packing fractions considered were in the range ζ ) 0.0001-0.53. The standard hard sphere MD procedure was used for the high-density systems, employing link lists to speed up the computations (e.g., ref 1, Chapter 3). At each density, typically ∼10 different values of N were simulated. The self-diffusion coefficient is the most straightforward of the hard sphere transport coefficients to compute because there is no “instantaneous” (i.e., t ) 0) contribution.18 Quantities quoted are in the usual hard sphere units; that is, the unit of mass is the mass of the sphere; m, the unit of energy

is kBT where kB is Boltzmann’s constant and T is the temperature; and the diameter of the hard sphere, σ, is the unit of length. The simulations were carried out with kBT ) 1, or T* ) 1 in more customary notation. At very low and decreasing packing fractions (e.g., , 0.1) the mean free path, λ, increases more rapidly than the side length of the simulation cell, as is now shown. From kinetic theory, λ/σ ) x2τ/Γ(1.5) where Γ(x) is the gamma function at x, and τ is the mean time between collisions for a particle, which is given by τ-1 ) 3(1 - N-1)(Z - 1)/xπ.17,7 For low packing fractions, Z - 1 = 4ζ, where Z is the compressibility factor.10 Therefore λ ∝ ζ-1 as ζ f 0, whereas the simulation cell sidelength L ) (πN/6ζ)1/3 ∝ ζ-1/3. This means that in the zero density limit (i.e., ζ f 0), the ratio of the mean free path to the simulation cell sidelength varies as λ/L f 0.2924ζ-2/3N-1/3, which, for example, is equal to unity when ζ ) 0.0071 for N ) 500. At very low packing fractions, it is, therefore, as expected, possible for the hard spheres to pass through the simulation cell without undergoing a collision. To accommodate this and to avoid missing collisions as a result of the periodic boundary conditions, the usual hard sphere MD procedure was modified.

D ) lim tc f ∞ of

1

N

∑ ∫0

3N i)1

tc

〈Vi(s + t)‚Vi(s)〉s dt

Hard-Sphere Fluid Self-Diffusion Coefficient

Figure 3. As for Figure 1, except that the diffusion data is plotted against N-2/3.

J. Phys. Chem. B, Vol. 111, No. 6, 2007 1457

Figure 5. The parameter A in eq 6 obtained by least-squares fitting the MD data to this formula for R ) 1/3. The symbols also have the vertical error bars on them. The line labeled “Theory 1/3” is eq 6 with the formula for A ≡ AT given below the equation, as derived by Fushiki14 and Yeh and Hummer.15 For AT, the MD shear viscosity, ηs, was used, as fitted to the formula ηs ) ηE f(ζ) where f(ζ) is given in eq 7 with the associated parameters provided below that equation. Note the data point with a value between 2 and 3 at ζ ) 0.01.

and was set to L/4 to make doubly sure. This procedure was found to be enacted only for packing fractions below 0.01, that is, at the very lowest of densities considered. All the collisions were found to be “real” for higher ζ. Note that this crossover packing fraction is at a much lower value than where the physical trends observed at low density (e.g., see Figure 7) first become evident. II. Theory and Simulation Results Kinetic theory for self-diffusion gives the self-diffusion coefficient in the limit of zero density,2 Figure 4. The ratio of the Enskog shear viscosity, ηE, from eq 17 in ref 10 to the MD generated shear viscosity, ηs, plotted against the packing fraction, ζ. The MD data was taken from ref 10, marked “MD” on the figure. The Erpenbeck and Wood N ) 4000 MD data from ref 25 are given in the Figure, as “EW” in the key. The MD data fitted to eq 7 are given in the Figure, and labeled “Fit”.

In the standard implementation of hard sphere molecular dynamics the scheme is said to be “event driven”, by which we mean that the program follows one collision after another in chronological order (rather than in time steps of fixed duration, as is used when the interaction potential is continuous). At very low densities, this procedure, unmodified, would lead to an incorrect collision sequence and even collisions being missed because the collision lists only take into account collisions with the nearest neighbors of a given particle. Because of the diverging mean free path, it is possible that the first collision a particle has with one of the other particles is with a nonnearest image. There are probably several ways of removing this problem. Here, a “virtual” or “null” collision was defined if the hard sphere was found to travel a certain distance before encountering another particle. The collision order lists were recomputed when either a real or virtual hard sphere collision had taken place first. By this procedure, the possibility of a nonnearest neighbor of a particle colliding first with the chosen particle was eliminated The distance a particle was allowed to travel before a virtual collision was deemed to have taken place has to be less than L/2 because of the nearest image convention

D0 )

( )

3σ kBT 8F πm

1/2

(3)

where F ) Nσ3/V, m is the mass of the hard sphere, kB is Boltzmann’s constant, and T is the absolute temperature. Note that D0 is inversely proportional to the density. DE is Enskog’s prediction for the self-diffusion coefficient.2

DE/D0 ) 1.01896/g(σ)

(4)

The radial distribution function at contact of the spheres, g(σ), can be written in terms of the compressibility factor, Z, g(σ) ) (Z - 1)/4ζ (e.g., see ref 10). There are many formulas for Z of hard spheres in the literature, and we use the analytic expression of Kolafa19,20 for Z

Z(ζ) )

1 + ζ + ζ2 - 2(ζ3 + ζ4)/3 (1 - ζ)3

(5)

because it is more accurate than the Carnahan-Starling formula at high density.21 At high packing fractions, deviations from the Enskog formula become apparent because of the importance of correlated collision sequences, which is why we have to resort to MD simulation to compute the D in this region. For thermodynamic quantities, all finite-size corrections to leading order vary as O(N-1) (e.g., see ref 22), as far as we are aware. The evidence to date from simulation studies for the self-diffusion coefficient (at not too low a density; see below)

1458 J. Phys. Chem. B, Vol. 111, No. 6, 2007

Heyes et al.

Figure 6. As for Figure 5 except that the low-density behavior of A is explored for R ) 1/3 and 1. The parameter A in eq 6 was obtained by least-squares fitting the MD data to this formula in the density region between ζ ) 10-4 and 0.1. The limiting density dependence of A(ζ) ) a1/ζa2, where a1 ) 0.0151 ( 0.0020 and a2 ) 1.080 ( 0.014 for R ) 1/3 and a1 ) 0.0959 ( 0.0136 and a2 ) 1.093 ( 0.015 for R ) 1. These “Fit” lines are given on the figure together with the MD values.

is for a significant but somewhat weaker system-size dependence. In this study, the set of MD self-diffusion coefficient data points computed for different Ns for a given ζ were fitted to the analytic form

D(ζ, N) ) D∞(ζ) - A(ζ) N-R

(6)

where D∞(ζ) is the value of the self-diffusion coefficient in the thermodynamic limit, and A(ζ) determines the magnitude of the finite size correction and its dependence on N. The formula given in eq 6 can be derived on the basis of hydrodynamics, in which case, R ) 1/3, and the analytic formula for the parameter A(ζ) in eq 6 is then A ) ξkBTζ1/3/62/3π4/3ηs, where the constant ξ = 2.837 297.15 ηs is the shear viscosity of the hard sphere fluid at the given packing fraction and temperature. Because kBT is a factorable scaling factor for all the transport coefficients of hard spheres, their temperature dependence is predictable without recourse to simulation investigation and, hence, can be omitted from this discussion. For future reference, the above analytic formula for A is denoted by AT. Here, R is generally treated as a free parameter in eq 6, in which case, formulas obtained for A(ζ) are empirical. Equation 6 has also been used to extrapolate to infinite N the self-diffusion coefficient of liquid para-hydrogen computed by ring-polymer molecular dynamics.23 In Figure 1, a selection of the D(ζ, N) values for the higher range of packing fractions is plotted against N-1/3, the procedure adopted in refs 14 and 15. They are fitted to the simulation data at each packing fraction using an implementation of the nonlinear least-squares Marquardt-Levenberg algorithm in the gnuplot graphics program.24 All data points are weighted equally because it is not easy to estimate the true errors in D for a given packing fraction and value of N. The apparent errors from the integral of the velocity autocorrelation function are usually quite small, certainly less than the symbol sizes on the figures. The main source of error, particularly at high fraction, probably arises from the limited sampling of phase space, which is difficult to quantify without performing even longer simulations. As may be seen, the data appear to conform to the N-1/3 scaling over the whole density range. Although the 1/3 exponent is expected on hydrodynamic grounds, this approximation only includes one

Figure 7. The least-squares best fit exponents R from the MD data using the formula in eq 6. The symbols have the vertical error bars with them. The curve on the figure is the best fit formula given in eq 8 applied to the data points on the figure.

possible source of the N dependence. One may question the general validity of the one-third exponent because other sources of N dependence may also be important, especially in the range of N accessible by simulation. Therefore, it is appropriate to investigate the quality of fit of the data to a more general plot of the form N-R, as in eq 6, where R is treated as an adjustable parameter at each density so that the function R(ζ) is obtained. Figures 2 and 3 show plots of this data with R ) 1/2 and 2/3, respectively. These figures show that R ) 1/2 appears to give as good a fit as R ) 1/3 for all N and packing fractions. The value R ) 2/3 is also quite a good fit at high packing fractions, but exhibits a noticeable upward pointing “cusp” for large N (on the left of the figure) for the lower ζ. An important point, therefore, is that the apparent quality of a fit guided by visual inspection is quite insensitive to not too extreme values of R. For the three values of R, both D∞ and A(ζ) increase with decreasing density. The importance of Figures 1-3 is that visual inspection with a preconceived value of the exponent is not reliable and that a nonlinear least-squares fit with the exponent treated as a variable is necessary. Since R is to be freely determined without bias, we must leave both it and A as free variables, both as a function of ζ. We can conveniently write the shear viscosity of the hard sphere fluid in the form ηs ) ηE f(ζ) where ηE is the analytic Enskog prediction for the shear viscosity (e.g., see ref 10), and the correction factor f(ζ) is established fitting to the MD data. f(ζ) f 1 for ζ f 0 as the Enskog prediction becomes increasingly more accurate. There have been a number of analytic expressions proposed in the literature for f but they are all functions that are valid only over certain packing fraction ranges, which makes them inconvenient for general analytical use (see, for example, the discussion in ref 10). Noting that the N dependence of the shear viscosity is much less than that of the self-diffusion coefficient, we find the MD shear viscosity data of Sigurgeirsson and Heyes can be fitted well within the statistical uncertainty for ζ < 0.54 with the new function,

f(ζ) ) 1.0 + a1ζa2

(7)

where a1 ) 2588 ( 258 and a2 ) 11.11 ( 0.13. The quality of fit can be seen in Figure 4, where ηE/ηs is plotted versus ζ. The Enskog formula is essentially exact for packing fractions less than ∼0.3. Also shown on the figure are the N ) 4000 MD

Hard-Sphere Fluid Self-Diffusion Coefficient

J. Phys. Chem. B, Vol. 111, No. 6, 2007 1459

TABLE 1: The Self-Diffusion Coefficient of the Hard Sphere Fluid in the Thermodynamic Limit, D∞, Obtained by Nonlinear Least-Squares Fit to the D(N) Simulation Data Using Eq 6, Where A and r Are also Treated as Free Parametersa ζ

D∞

σD

ζ

D∞

σD

ζ

D∞

σD

0.0001 0.001 0.005 0.010 0.050 0.100 0.150 0.200 0.250

1127 112.5 22.4 11.05 2.039 0.936 0.5785 0.402 0.287

3 0.6 0.2 0.06 0.006 0.001 0.004 0.004 0.003

0.275 0.300 0.325 0.350 0.375 0.400 0.410 0.420 0.430

0.240 0.2021 0.168 0.1396 0.1136 0.0877 0.0790 0.0699 0.0620

0.002 0.0009 0.003 0.0012 0.0066 0.0006 0.0009 0.0007 0.0008

0.440 0.450 0.460 0.470 0.480 0.490 0.500 0.510

0.0549 0.0485 0.0413 0.0351 0.0301 0.0249 0.0201 0.0153

0.0007 0.0003 0.0005 0.0005 0.0004 0.0002 0.0002 0.0002

a

The units of D∞ are σ(m/kBT)1/2; σD is the standard error.

data from Erpenbeck and Wood,25 which agrees within statistics with the MD data given in ref 10. Figure 5 shows the packing fraction dependence of the coefficient A in eq 6 with R ) 1/3. The A’s are obtained by nonlinear least-squares fit to the MD data. The analytic function AT is also shown in the Figure. For each density in Figure 1, the parametric errors in the values of A from the nonlinear leastsquares fit were obtained, and these are shown in Figure 5. Agreement between theory and simulation data is good in the approximate packing fraction range 0.15-0.4; however, for lower packing fractions, AT decreases with diminishing ζ, whereas the MD data indicates an increase in this quantity. This trend is independent of the value of R chosen. In fact, the MD data suggest a limiting form ∼A ∝ ζ-1 for low densities (see Figure 6, which focuses on the low density behavior of A), just as for the mean free path. The AT function also appears to underestimate A for packing fractions in excess of ∼0.4. We suggest that in the two density limits, there are important additional sources of N dependence in D in the covered range of N. At low packing fraction, the mean free path (∼ζ-1) increases more rapidly with decreasing density than the simulation side length (∼ζ-1/3) so that the effects of the periodic boundaries on D for a given value of N increase as density decreases. In fact, Alexander et al.26 have already noted using direct simulation Monte Carlo, DSMC, that significant changes of the transport coefficients occur when the cell dimensions are larger than the mean free path. DSMC is a stochastic simulation technique that solves the Boltzmann equation by replacing the distribution function with a representative set of particles.27 At high packing fractions, excluded volume effects become significant. In addition, another potential source of an N dependence of D was discovered by Alder and Wainwright.3,28 A study of the velocity autocorrelation function for hard spheres over the entire fluid range showed a deviation from the exponential relaxation expected from Enskog theory. An algebraic decay with time of the VACF whose effects were most prominent at intermediate packing fractions = 0.24 was observed. The resulting diffusion coefficient was up to ∼16% larger than the Enskog prediction in this region. They showed that this was due to an enhanced cooperative movement of the trajectories of the hard spheres, whose average effects can be accounted for using the hydrodynamic continuity equations.28 This important local hydrodynamic contribution to the selfdiffusion coefficient had not been predicted before that simulation and theoretical study. Therefore, the system size dependence of the self-diffusion coefficient across the fluid density range is more complicated than the hydrodynamic formulas of refs 13 and 14 might lead one to expect. Nevertheless, the formula of Dunweg and Kremer13 is significant in being the first entirely analytic prediction of D(N) that accounts very well for the simulation data at not too high or too low densities. In addition, as Yeh and Hummer15 demonstrated, their formula is not

restricted to hard spheres and should be applicable to any simple molecular liquid. Figure 7 shows the density dependence of the exponent, R in the formula of eq 6. Despite the significant amount of selfdiffusion coefficient data, the values of R are quite scattered. Nevertheless, a trend in R can be discerned in which it apparently decreases to ∼1 as ζ f 0, goes through a minimum of R = 1/3 at ζ = 1/3, and then increases to = 2/3 or possibly even higher at ζ = 0.5. A polynomial fit to the data and shown in Figure 7 is

R(ζ) ) 1 + a1ζ + a2ζ2 + a3ζ3 + a4ζ4 + a5ζ5

(8)

where a1 ) -4.493 57, a2 ) 9.215 73, a3 ) 4.520 95, a4 ) -49.6403, and a5 ) 69.0057. Note that this assumes R ) 1 in the zero density limit. The data error estimates were used to calculate the relative weight of each data point when determining the weighted sum of squared residuals for the nonlinear leastsquares fit. The standard deviation of the ordinate value, σ, is used to compute a weight for the datum, σ - 2. Substitution of eq 8 in eq 6 gives, on least-squares fitting of the MD self-diffusion coefficients, an optimized set of D∞(ζ) and A(ζ) determined with the most likely value of R at each packing fraction. An even simpler analytic form, which is almost as good a fit to the simulation exponents, is to choose R(ζ) ) 1 + e1ζ + e2ζ2 and assume the function goes through the (R, ζ) points (1, 0), (1/3, 1/3), and (2/3, 1/2) which gives e1 ) -14/3 and e2 ) 8. The analytic theory leading to the hydrodynamic or N-1/3 scaling is probably only asymptotically correct for large N close to the thermodynamic limit. Other mechanisms give rise to a more complicated dependence of the self-diffusion coefficient on N and packing fraction. These could be accommodated in the more general expression DN(N, ζ) ) ∑R cR(N, ζ)/NR, for example, where the functions cR(N, ζ) are used to associate the different mechanisms to certain regions of N and ζ. The R ) 1/3 value would be the leading term and dominant for large N and not too low density. Larger R values are important at high packing fractions and small N values, which derive from excluded volume effects. For low densities, particularly large N values were not used in the simulations; for example, N ) 864 was the maximum for a packing fraction of 0.001, and for ζ > 0.4, much larger N (up to 10 976) were used. If the fitting is restricted to N < 300 in the high-density region, the R values are larger than if all the N values are taken into account (as is the case for Figure 7). The value of R increased to ∼1 at a ζ ) 0.5 for N < 300. The reason the fitted exponent function shown in Figure 7, obtained from all N simulation data, was used to obtain D∞ and A in the thermodynamic limit is that the largest N values used in the simulations are not obviously in the R ) 1/3 scaling regime. Even if they were, there would not be enough points to achieve satisfactory statistics using the R ) 1/3 value.

1460 J. Phys. Chem. B, Vol. 111, No. 6, 2007

Heyes et al.

Figure 8. The ratio D∞/DE from the MD simulations based on eq 6 using the analytic formula for R(ζ) given in eq 8. The errors in the data points from the extrapolation are typically less than ∼1% and within the scale of the symbol. The continuous curve (“Fit”) is the formula for D∞/DE taken from eq 11. The curve “EW” is the formula of Erpenbeck and Wood for the self-diffusion coefficient in the thermodynamic limit given in eq 10.

The extrapolation to the thermodynamic limit used here is based on an “effective” or average value of R applied to the whole N range used at each packing fraction. This should be a more accurate route to the thermodynamic limit because there are larger statistical errors associated with just using the few largest N points for extapolation to infinite N, even though these points are closer to the thermodynamic limit. Table 1 and Figure 8 show the extrapolated values of D∞ resulting from this procedure as a function of ζ. The ratio D∞/ DE is shown in the figure, which is unity in the zero density limit, and it increases with density. A value of D∞/DE in excess of unity indicates an enhanced self-diffusion coefficient above that of the Enskog prediction, a phenomenon first discovered by Alder and Wainwight in their seminal 1967 publication.3 For packing fractions in excess of ζ ) 0.30, the ratio D∞/DE decreases and falls below 1 for packing fractions greater than ∼0.43, which reflects “cage” entrapment of the particles. These trends were encapsulated for the first time in a simple formula by Speedy who used the hard sphere MD simulation selfdiffusion data available at the time,8

F ( (1.09 ))(1 + F (0.4 - 0.83F ))

D/D0 ) 1 -

2

2

(9)

where D0 is the self-diffusion coefficient in the limit of zero density, defined in eq 3. The first factor in eq 9 represents the cage constraining effects at high density, which leads eventually to a glass transition (i.e., where D f 0) at F ) 1.09 or ζ ) 0.571. The second term accounts for the hydrodynamic enhancement of D associated with an algebraic decay of the VACF (see below), discovered by Alder and Wainwight. Erpenbeck and Wood6 fitted their MD hard-sphere, self-diffusion coefficient simulation data when extrapolated to infinite N by the following formula, normalizing this time by the Enskog self-diffusion coefficient,

() () ()

D∞/DE ) 1 + a1

V0 V0 2 V0 + a2 + a3 V V V

3

(10)

where DE is Enskog’s prediction for the self-diffusion coefficient,2 defined in eq 4. In the above formula, V is the volume available to each particle, and a1 ) 0.054 034 49, a2 )

Figure 9. As for Figure 8, except that -Aζ is plotted against ζ. The symbols also indicate the error bars in ordinate values. The continuous curve is taken from eq 12.

6.365 616, and a3 ) -10.942 539. The quantity V0 is the closepacked volume per particle for a fcc crystal which is equal to σ3/x2, and V0/V ) 6ζ/πx2. The formula in eq 10 gives D∞ ) 0 at ζ ) 0.5557. The Enskog normalization is also preferred here because it accounts for the MD data to a higher packing fraction than simply using D0 in the denominator. The selfdiffusion coefficients in the thermodynamic limit, that is, N f ∞, D∞, computed in the present study were fitted by nonlinear least squares to the following analytic form,

D∞/DE ) 1 + b1ζ + b2ζ2 + b3ζ3 + b4ζ4 + b5ζ5 (11) where b1 ) 0.277 645, b2 ) 3.989 64, b3 ) 26.498 81, b4 ) -134.0015, and b5 ) 110.7344. The D∞ were obtained from eq 6 with R calculated using eq 8. The formula in eq 11 gives D∞ ) 0 at ζ ) 0.5602. Figure 8 compares the expression in eq 11 with the Erpenbeck and Wood (EW) formula of eq 10.6 The agreement is generally quite good, although the present formula predicts a slightly lower value for D∞ than that of EW for densities in excess of ∼0.40. In fact, statistically, the same values between the two formulas for ζ > 0.40 can be obtained with a fixed value of R ) 1/3 at all packing fractions (in which case, b1 ) 1.586 46, b2 ) -3.505 39, b3 ) 30.5922, b4 ) -94.4749, and b5 ) 55.3039 in eq 11). As we have discussed, eq 6 is a convenient analytic form with which to fit the data, but it is probably not unique. Figure 9 shows the simulation-derived A(ζ), again using the R(ζ) function from eq 8 substituted in eq 6. A least-squares fit to these values is given in the Figure, which has the form,

-Aζ ) c1 + c2ζ + c3ζ2 + c4ζ3 + c5ζ4

(12)

and where c1 ) 0.183 724, c2 ) -1.304 907, c3 ) 6.309 320, c4 ) -15.370 63, and c5 ) 13.701 86. The function A diverges as ζ f 0, and decreases with increasing ζ, with some evidence of an increase again at the very highest densities. The trend in R is, therefore, mirrored by that in A. To summarize, the packing fraction and N-dependent, selfdiffusion coefficient, D(ζ, N), can be accounted for within statistical error using eq 6, with D∞ given by eq 11, A(ζ) given by eq 12, and R(ζ) from the expression in eq 8. D(ζ, N) curves for selected values of N are given in Figure 10 and compared with values taken from the literature from other MD studies. The MD data from Lue17 for N ) 1000 is reproduced well by this parametrization of eq 6. Note that the maximum sensitivity

Hard-Sphere Fluid Self-Diffusion Coefficient

J. Phys. Chem. B, Vol. 111, No. 6, 2007 1461

Figure 10. The MD-generated self-diffusion coefficients compared with the Enskog value as a function of ζ for a series of N values. The MD N ) 108 data points (“108” on the figure) are those given in Figure 1. The continuous curves are calculated for the various N values using eq 6 with D∞ given by eq 11; A(ζ) by eq 12; and R(ζ), by eq 8. The data “Lue 1000” is taken from the N ) 1000 MD study of ref 17, data given in Table 1 of that paper. “EW 4000” and “EW N ) ∞” are the MD data and extrapolations from the Erpenbeck and Wood MD study.6

to N is for a packing fraction of ∼0.3 and that the self-diffusion coefficient for N ) 108, say, is about 15% lower than the estimated value at the thermodynamic limit at this packing fraction. Also shown in the Figure are the MD data for 4000 hard spheres from Erpenbeck and Wood6 and their extrapolation to N f ∞. These data agree quite well with the corresponding predictions of the parametrized formula in eq 6. Using eq 11 for the diffusion coefficient (i.e., in the thermodynamic limit), there is a maximum in the ratio at a packing fraction of ζmax ) 0.2971 with a value of D∞/DE ) 1.346 94, which compares favorably with the values ζmax ) 0.2945 and D∞/DE ) 1.342 from the formula in eq 10 of Erpenbeck and Wood.6 The formula of Speedy from eq 9 gives the somewhat smaller values of 0.2809 and 1.232, respectively, which one can attribute to the relatively small values of N used for his fit. If we employ our prescription for D(ζ, N), using eq 6, with D∞ given in eq 11, A(ζ) given in eq 12 and R(ζ) in eq 8, then for N ) 108, 256, 500, 864, and 10 976, there is ζmax ) 0.2647, 0.2728, 0.2773, 0.2803, and 0.2884, and the corresponding maximum values of D∞/DE are 1.131, 1.188, 1.222, 1.244, and 1.303, respectively. The Speedy formula agrees best with the present data for N = 1000. An alternative parametrization of D∞/DE from eq 11 is to use an expression similar to that of Speedy.8

(

D∞/DE ) 1 -

)

ζ d1 (1 + d2ζ + d3ζ2 + d4ζ3) ζg

(13)

where a nonlinear least-squares fit gives, ζg ) 0.5531, and d1 ) 0.041 280, d2 ) -0.160 226, d3 ) 13.8375, and d4 ) -30.3387. Note that the self-diffusion coefficient in the thermodynamic limit in this formula is normalized by the Enskog prediction for the self-diffusion coefficient rather than by D0, as in ref 8. The advantage of the analytic form given in eq 13 is that it attempts to separate the crowding effect (first factor) from the Alder-Wainwright local scale hydrodynamic effect (second factor) in eq 13. This formula reproduces the D∞ in Figure 8 within the statistical uncertainty of the simulation data.

Figure 11. The velocity autocorrelation function of the hard sphere fluid at a packing fraction of 0.3. The curve labeled “Enskog” is taken from eq 17, “Long tail” from eq 18 multiplied by a factor of 0.4. The other curves are the simulation data for different N values.

The velocity autocorrelation function from Enskog theory of hard spheres is an exponential,2

( )

CV(t) ) exp -

2t 3τE

(14)

where

τE )

xπσ 6(kBT/m)1/2(Z - 1)

(15)

Therefore, on substituting eq 15 into eq 14,

( ( )( ) )

CV(t) ) exp -4

Z - 1 kBT σ πm

1/2

t

(16)

where the compressibility factor can be calculated from, for example, eq 5. Figure 11 shows the VACF at ζ ) 0.3 for N values between 32 and 10 976. As the system size increases, the VACF decays more slowly at long times and develops a linear behavior on the log-log plot indicative of an algebraic with time decay of the VACF. The various contributions to the velocity autocorrelation function of the hard sphere fluid arising from correlated particle motion were discussed by Erpenbeck and Wood6 (see also refs 29 and 30). The Enskog prediction for the VACF given in eq 16 is excellent at short time for all densities. At ζ ) 0.3, the two curves are statistically indistinguishable down to a value of ∼0.6 at a time of 0.075. Differences between these two curves arising from correlated particle motion become more important with time. At long times, the dominant contribution in this category comes from the shear modes, which gives a long time algebraic tail to the VACF (see also ref 31),

CV(t) )

RD t3/2

(17)

where

RD )

2kBT 3Fm[4π(D + ν)]3/2

(18)

1462 J. Phys. Chem. B, Vol. 111, No. 6, 2007

Heyes et al.

Figure 12. The boundary condition parameter, c ) kBT/πD∞ηs, as a function of packing fraction. D∞ was obtained from eq 11, and the shear viscosity, from ηs ) ηE f(ζ), where f(ζ) is given in eq 7. Note that c ) 2 corresponds to a slip boundary condition, whereas c ) 3 is for the “stick” limit.

F ) N/V and ν ) ηs/mF is the kinematic viscosity. Therefore, assuming this form for CV for t > tc, the long time correction to the self-diffusion coefficient is 2RDtc-1/2. The value of tc was taken to be at least several hundred collisions per particle and for the highest densities, tc = 100 hard sphere time units. The long time tail corrections are from 2RDtc-1/2 less than ∼1% (for the highest density, taking into account that the theoretical value for RD overestimates the computed value from eq 18), and its relative contribution decreases dramatically with decreasing density. In addition, as was mentioned earlier, the VACFs were observed in all cases to decrease statistically to zero before truncation. This local hydrodynamic enhancement to D appears to be significant only at intermediate densities. In the very lowdensity and high-density limits, its effects are insignificant. There are in the literature various semiempirical relationships between the self-diffusion coefficient and the shear viscosity, notably the Stokes-Einstein, SE, formula,

Dηs )

kBT cπσ

(19)

where the boundary condition parameter, c, is 3 and 2 for “stick” and “slip” boundary conditions, respectively. Although SE strictly only applies to macroscopic spheres in a Newtonian liquid, it has been tested for dense liquids and found to work remarkably well on the molecular scale with some flexibility in assigment of the molecular diameter or c parameter (see, e.g., ref 32, p 262). The optimum value of c is found to be between the stick and slip limits, depending somewhat on the state point.33 For molecular fluids, there are variants on eq 19 in the literature, in which the fluid density is included; for example, see refs 34-36. This is reasonable because in the zero density limit, ζDηsσ/kBT f 5/256 from kinetic theory, using definitions for the self-diffusion coefficient and shear viscosity of hard spheres in the zero density limit, given in, for example, ref 10. However, for hard spheres using the present data, we find that the quantity, 256ζD∞ηs/5 increases monotonically from unity in the zero density limit to a value of ∼3.8 at ζ ) 0.5. The hard sphere data of this work, extrapolated to the thermodynamic limit, reveals a nontrivial density dependence of D∞ηs, which can be represented in terms of an effective boundary condition coefficient, c(ζ), as shown in Figure 12.

Figure 13. The product Dηsp against the packing fraction, where p ) 0.975. The figure shows that the relationship, Dηs0.975 ) 0.15 applies well over most of the equilibrium fluid regime, that is, for packing fractions between ∼0.2 and 0.5.

(For a recent discussion of the application of the SE relationship on the molecular scale see, e.g., ref 37.) This figure shows that c increases to a value just above 2 for ζ > 0.2, with a further dramatic increase for ζ > 0.5. One cannot be sure at this stage whether the ζ > 0.5 behavior shown here is a “real” effect because the fluid is in a metastable state in this packing fraction regime, and a more extensive simulation study in this regime is required. Indeed, the fluctuation-dissipation formula for the self-diffusion coefficient (i.e., eq 2) is not strictly valid in a nonequilibrium state, such as a metastable fluid, that is, for ζ > 0.493 here9 (see e.g., refs 38 and 39 for a discussion of violation of the fluctuation-dissipation relation). An alternative proposal is that Dηsp ) k, where k is a constant, and p = 0.686.40-42 Treating the exponent p as an adjustable parameter, Figure 13 demonstrates that this relationship, with p ) 0.975, is constant within the computational uncertainties for 0.2 < ζ < 0.5, and we find k = 0.15 in this packing fraction range. Another approach to predict the self-diffusion coefficient is to express it in terms of the excess entropy per particle, sex, where43

sex/kB ) -

∫0ζ Z -ζ 1 dζ

(20)

where Z is the compressibility factor, which can be evaluated analytically for a given equation of state analytic in ζ. The idea that entropy can be used to scale the transport coefficients lies in its being a measure of the confinement of molecules by their neighbors, which should affect both sets of properties. For the Carnahan-Starling21 and Kolafa19 equations of state,

sex/kB ) -ζ

(4 - 3ζ) (1 - ζ)2

(21)

and

(4ζ2 - 33ζ + 34) 5 sex/kB ) - ln(1 - x) - ζ 3 6(1 - ζ)2

(22)

respectively (the latter being considered new). A number of semiempirical formulas where D ) a exp(bsex/kB), in which the parameter a is usually density- and temperature-dependent and b is a constant (e.g., see refs 44-47), have been proposed

Hard-Sphere Fluid Self-Diffusion Coefficient

Figure 14. A plot of ln(D∞/DE g(σ)2)/γ, where γ ) 0.8, from the extrapolated MD data, labeled “MD N ) ∞” on the figure. Also shown are sex/kB after Bastea.46 The excess entropies, sex, were computed from the formulas given in eqs 21 for the Carnahan-Starling equation of state and eq 22 for the Kolafa equation of state, respectively. On the Figure, it is hard to distinguish the latter two curves

in the literature for hard spheres. One of the more successful is D∞/DE ) g(σ)2 exp(γsex/kB) by Bastea,46 where g(σ) is the value of the radial distribution function at contact and γ ) 0.8, which is tested in Figure 14. (Note that the excess entropies computed using the formulas in eqs 21 and 22 are essentially indistinguishable on the plot.) The Bastea formula reproduces the simulation-derived D∞ very well up to a packing fraction of ∼0.45. It might be expected that such a treatment in terms of the entropy is best founded at high packing fractions (e.g., ∼ζ > 0.4), where particle confinement is the dominant factor in governing D. In the present context, perhaps the most significant aspect of any potential relationship between self-diffusion and excess entropy is that thermodynamic properties are known to have an N-1 scaling (e.g., see ref 22), which would imply the same for the self-diffusion coefficient. This is consistent with the increasing value of the exponent R toward 2/3 and possibly higher in the high-density regime, as revealed in Figure 7. III. Conclusions To conclude, molecular dynamics simulations have been carried out to determine the self-diffusion of the hard sphere fluid for a wide range of packing fractions and system sizes. The data have been fitted to a simple analytic form that is a generalization of one derived by Dunweg and Kremer,13 Fushiki,14 and Yeh and Hummer15 in eq 6, where D is linear in N-1/3. As a word of warning, just because a plot of D against N-R for a particular value of R gives an apparently good, straight line does not mean to say that the R value is optimum. Indeed, for the relatively small system sizes accessible by simulation, the apparent quality of fit to the eye is quite insensitive to fairly large variations in R. Even with modern computer resources, it is still problematic to characterize the system size dependence of the self-diffusion coefficient. In the present treatment for optimum effect, R was chosen to be an adjustable parameter, and a least-squares fit to the data revealed a strong dependence of its value on the packing fraction, ζ. For ζ = 0.35, it is found that R = 1/3, but is larger than 1/3 at both low and high densities. The N dependence is probably the result of a combination of effects whose relative importance depends on density and the range of N values studied. For low packing fractions, for ∼ζ < 0.1, the divergence

J. Phys. Chem. B, Vol. 111, No. 6, 2007 1463 in the mean free path is probably the most important factor, leading to R = 1 in the zero density limit, whereas at intermediate densities, macroscopic hydrodynamic effects are dominant (giving R = 1/3). The analytical theory is just based on hydrodynamics, which might explain why it does not fully account for the low-density data. At higher packing fractions in the studied range of system sizes, excluded volume factors probably dominate so that the D scaling follows more closely that already characterized for the thermodynamic properties (see, e.g., Kolafa et al. in ref 22). It is, of course, possible that for a given packing fraction, the proposed analytic form for D(ζ, N) is not as given in eq 6, but of a more general form, DN(ζ, N) ) ∑RcR(N, ζ)/NR; for example, where the functions cR(N, ζ) are used to associate the different mechanisms to certain regions of N and ζ, where the N dependence for small N values is especially strong at high densities. The curve obtained here for the self-diffusion coefficient in the thermodynamic limit, derived from a density dependent R, follows closely that obtained by Erpenbeck and Wood for ζ < 0.4, but the D values are a little lower at higher packing fractions than the EW prediction. Acknowledgment. The authors thank M. Fushuki for supplying tabulated data from the paper of ref 14 and the referees for helpful comments and suggestions. References and Notes (1) Hansen, J. P.; McDonald, I. R. Theory of Simple Liquids, 2nd ed.; Academic Press: London, U.K., 1986. (2) Chapman, S.; Cowling, T. G. The Mathematical Theory of Nonuniform Gases, 2nd ed.; Cambridge University Press: Cambridge, U.K. 1952. (3) Alder, B. J.; Wainwright, T. E. Phys. ReV. Lett. 1967, 18, 988. (4) Alder, B. J.; Gass, D. M.; Wainwright, T. E. J. Chem. Phys. 1970, 53, 3813. (5) Dymond, J. H.; Alder, B. J. J. Chem. Phys. 1968, 48, 343. (6) Erpenbeck, J. J.; Wood, W. W. Phys. ReV. A 1991, 43, 4254. (7) Woodcock, L. V. Ann. N. Y. Acad. Sci. 1981, 371, 274. (8) Speedy, R. J. Mol. Phys. 1987, 62, 509. (9) Speedy, R. J. Mol. Phys. 1998, 95, 168. (10) Sigurgeirsson, H.; Heyes, D. M. Mol. Phys. 2003, 101, 469. (11) Liu, H.; Silva, C. M.; Macedo, E. A. Chem. Eng. Sci. 1998, 53, 2403. (12) Kob, W.; Andersen, H. C. Phys. ReV. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top. 1993, 48, 4364. (13) Dunweg, B.; Kremer, K. J. Chem. Phys. 1993, 99, 6983. (14) Fushiki, M. Phys. ReV. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top. 2003, 68, 021203. (15) Yeh, I. C.; Hummer, G. J. Phys. Chem. B 2004, 108, 15873. (16) Heyes, D. M. The Liquid State; John Wiley & Sons: Chichester, U.K., 1997. (17) Lue, L. J. Chem. Phys. 2005, 122, 044513. (18) Ernst, M. H.; Brito, R. Phys. ReV. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top. 2005, 72, 061102. (19) Boublik, T.; Nezbeda, I. Coll. Czech. Chem. Commun. 1986, 51, 2301. (20) Ayala de Lonngi, D.; Lonngi Villanueva, P. A. Mol. Phys. 1991, 73, 763. (21) Carnahan, N. F.; Starling, K. E. J. Chem. Phys. 1969, 51, 635. (22) Kolafa, J.; Labik, S.; Malijevsky´, A. Phys. Chem. Chem. Phys. 2005, 6, 2335. (23) Miller, T. F., III; Manolopoulos, D. E. J. Chem. Phys. 2005, 122, 184503. (24) Williams, T; Kelley, C. GnuplotsAn interactive plotting program, Version 3.7; organised by D. Denholm, 1998. (25) Erpenbeck, J. J.; Wood, W. W. Phys. ReV. A 1985, 32, 412. (26) Alexander, F. J.; Garcia, A. L.; Alder, B. J. Phys. Fluids 1998, 10, 1540. (27) Bird, G. A. Molecular Gas Dynamics and the Direct Simulation of Gas Flows; Clarendon Press: Oxford, U.K., 1994. (28) Alder, B. J.; Wainwright, T. E. Phys. ReV. A: At., Mol., Opt. Phys. 1970, 1, 18. (29) McDonough, A.; Russo, S. P.; Snook, I. K. Phys. ReV. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top. 2001, 63, 026109.

1464 J. Phys. Chem. B, Vol. 111, No. 6, 2007 (30) Dorfman, J. R.; Kirkpatrick, T. R.; Sengers, J. V. Ann. ReV. Phys. Chem. 1994, 45, 213. (31) Widom, A. Phys. ReV. A: At., Mol., Opt. Phys. 1971, 3, 1394. (32) Balucani, U.; Zoppi, M. Dynamics of the Liquid State; Clarendon Press: Oxford, U.K., 1994; p 18. (33) Heyes, D. M.; Bran´ka, A. C. Phys. Chem. Chem. Phys. 2005, 7, 1220. (34) March, N. H.; Tosi, M. P. Phys. ReV. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top. 1999, 60, 2402. (35) Zwanzig, R. J. Chem. Phys. 1983, 79, 4507. (36) Rah, R; Eu, B. C. Phys. ReV. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top. 1999, 4, 4105. (37) Stillinger, F. H.; Debenedetti, P. G. J. Phys. Chem. B 2005, 109, 6604.

Heyes et al. (38) Berthier, L.; Barrat, J. L. J. Chem. Phys. 2002, 116, 6228. (39) Bonn, D.; Kegel, W. K. J. Chem. Phys. 2003, 118, 2005. (40) Zwanzig, R.; Harrison, A. K. J. Chem. Phys. 1985, 83, 5861. (41) Chauhan, A. S.; Ravi, R.; Chhabra, R. P. Chem. Phys. 2000, 252, 227. (42) Ould-Kaddour, F.; Barrat, J. L. Phys. ReV. A: At., Mol., Opt. Phys. 1992, 45 2308. (43) Senthil Kumar, V.; Kumaran, V. J. Chem. Phys. 2005, 123, 114501. (44) Rosenfeld, Y. J. Phys. Condens. Matt. 1999, 11, 5415. (45) Bretonnet, J. L. J. Chem. Phys. 2002, 117, 9370. (46) Bastea, S. Phys. ReV. Lett. 2004, 93, 199603. (47) Samanta, A.; Musharaf Ali, Sk.; Ghosh, S. K. Phys. ReV. Lett. 2004, 93, 199604.