Rotational Relaxation in Liquids Molecular Hydrodynamics and Long Time Tails BY BRUCEJ. BERNE Department of Chemistry Columbia University New York N.Y. 10027 U.S.A. Received 24th September 1976 1. MOLECULAR HYDRODYNAMICS It is remarkable that hydrodynamic models can be applied with great accuracy at the molecular level. The evidence for this comes from several quarters. For example Alder et d.,'using computer experiments have shown that the translational friction constant c of a sphere in a one-component dense smooth hard sphere fluid is c = 4nqa where 11 is the shear viscosity and a is the radius of the sphere. This result is identical to Stoke's calculation of c for asphere with slip boundary conditions. Recently a variation principle has been used to evaluate the rotational friction co- efficient of a uniformly rotating slippery spheroid and has shown that this can differ from the friction coefficient of a sticky spheroid by as much as an order of magnitude or more.Assuming that benzene is a slippery ellipsoid Zwanzig and Hu's calculation is in excellent agreement with rotational diffusion coefficients determined by de- polarized light ~cattering.~ It would appear from the foregoing that the details of the intermolecular potential are unimportant when it comes to the determination of translational and rotational friction constants except insofar as these details effect the value of the viscosity. The success of hydrodynamics does not stop with the determination of static friction constants. Several years ago Zwanzig and Bixon calculated the normalized velocity correlation function by solving the Navier Stokes equation in the low Reynolds number limit for a particle executing non-uniform translatory motion in a visco-elastic continuum fluid with boundary conditions intermediate between pure stick and pure slip boundary condi- tions.This model has the general characteristic that at long times the correlation function decays asymptotically as a,t-3/2 where the coefficient a and the exponent are in excellent agreement with the asymptotic time dependence observed by Alder and Wainwright in their computer experiments on smooth hard sphere fluids. More- over Verlet et a1.' have shown that the hydrodynamic model gives an excellent fit to the detailed structure of the velocity correlation function determined by them for a dense Lennard-Jones fluid.Several years ago Ailawadi and Berne using hydrodynamic equations containing the angular momentum variables showed that the normalized angular velocity correla- tion function BRUCE J. BERNE decays asymptotically as C"(t) = badt-'d +2)/2 (103) where and where I my n v D,and d are respectively the molecular moment of inertia molecular mass number density kinematic shear viscosity (v = q/mn) self diffusion coefficient and dimensionality of the system. This should be compared with the asymptotic prediction for the linear velocity cv(t) = avdf-d'2 (1 5) where -1 avd = 4nd -[4n(D + v)]"'2. The latter result has been confirmed by computer experiment for d = 2 3.Berne9 subsequently calculated the full time dependence of Cm(t)for a sticky sphere executing non-uniform rotations in a viscous continuum fluid and showed that in the limit t -+ 00 it reduces to eqn (1.3). Eqn (1.3) has also been confirmed by mode-mode calculations. Assuming the validity of a generalized rotational diffusion equation Berne conjectured that the orientational correlation functions G(t>=(Pl(U(0) u(t))> (1 -7) also behave asymptotically as This asymptotic behaviour also follows directly from the "hydrodynamic equa- tions " that have been used to calculate the so called " shear doublets "in depolarized light scattering and thus appears to be quite reasonable.lo Subsequent mode-mode calculations l1 have taken issue with eqn (1.7) and there exist several other predictions for the asymptotic form of Cl(t).Recently Hill and Deutch,12 using Faxen's theorem have shown that eqn (1.7) is valid. They conjecture that the velocity fields arising from finite sized particles give eqn (1.8) whereas the mode-mode coupling theory calculations correspond to velocity fields arising from point sources. Recently Montgomery and Berne13 have computed Cco(t)for a sticky sphere executing non-uniform rotations in a viscoelastic continuum fluid. This model gives the asymptotic result stated in eqn (1.3) no matter what model is used for the visco- elasticity. It is shown nevertheless that the short time behaviour is strongly depend- ent on the details of the viscoelastic model.For example (1.9a) L (1 -nl/r t + o(t2). viscoelastic (1.9b) It is interesting to compare this hydrodynamic model for Co(t)with the rough ROTATIONAL RELAXATION IN LIQUIDS sphere molecular dynamics of O’Dell and Berne. l4 First note that the initial decay of Co(t)in a rough sphere fluid is given by (1.10) where K is the dimensionless moment of inertia of a sphere of diameter 0 (K = 41/ma2) and tEis the Enskog relaxation time. Comparison of eqn (1.10) and (1.9) shows that if the hydrodynamic model is to describe the rough sphere fluid a viscoelastic theory must be used instead of the viscous theory. Eqn (1.10) can then be used to eliminate one of the two parameters in the model. Fig 1.1 gives a comparison between the “best O00,2 4 6 .8 10 12 *1 ~ time in mean collision times FIG.1.1 .-The dots ...indicate the values of Cw(t)determined from computer experiments and the solid line indicates the “ best fit ” theory.fit ” single relaxation time viscoelastic model and a computer experiment on a rough sphere fluid containing 108 spheres and undergoing 250 OOO collisions. The results look very encouraging and we might conclude that the hydrodynamic model is excellent. Nevertheless the comparison involves the choice of a single parameter. A further check can be made. We note from eqn (1.3) that the model gives a long time tail which depends on D and v thus we can relate the parameter to be determined to v and D. We find that the best fit parameter differs dramatically from the value calculated for it using the Enskog values of v and D.Thus there is no internal consistency. Putting it another way if we determine the parameter from the Enskog values of v and D,we can determine Co(t)on the basis of our hydro- dynamic model. Fig. 1.2 shows a comparison between this procedure [curve @)] and experiment [curve (a)]. The troublesome thing about the hydrodynamic theory is that if it is forced to give the supposedly correct initial and long time decay it predicts strong oscillations in Co(t). Although oscillations like this sometimes appear in fluids containing structured particles they do not appear to occur in rough sphere fluids. Because the spheres are microscopically rough we made the plausible assumption that stick boundary conditions should apply.In hydrodynamics the rotating particle creates a boundary layer. This boundary layer can be accounted for by introducing a BRUCE J. BERNE FIG.1.2.-Curve (a)is a result of the computer experiment and curve (b)is C,(t) as determined from the model by fitting the initial decay to eqn (1.10)and long time decay to eqn (1.3) with Y,D calculated from Enskog theory. parameter of slip p. To allow for arbitrary slip at the surface of the rotating sphere the boundary condition is generalized so that the tangential component of the force is proportional to the velocity of the fluid relative to the velocity of the particle surface. a,. a.a = B [u -S2 x R]. a,. Here p = 00 corresponds to the complete stick limit discussed above and /? = 0 corresponds to pure slip (i.e.free rotation). This is discussed in ref. (13b) with ap- propriate references. A straightforward calculation gives the following short time behaviour 1-2Pt + O(P) viscous (1.lla) ~w(0 tz t + o(t2)viscoelastic one relaxation time. (1.1 1b) Comparison of eqn (1.1 la) with eqn (1.9~)shows that slip boundary conditions drastic- ally alter the short time behaviour so that the viscous model is now consistent with the exact short time decay expected from kinetic theory. Since L is known and q is given by kinetic theory it is possible to compute p by forcing eqn (1.1 la) to agree with eqn (1 .lo). On the other hand since we do not know y we cannot evaluate B from eqn (l.llb)andeqn(l.lO). An alternative viscoelastic model is given by Theodosoplu and Dahler from the kinetic theory of a dense rough sphere fluid with q yl and y2 expressed in ROTATIONAL RELAXATION IN LIQUIDS terms of microscopic quantities.The short time dependence of Co(t)for this model is identical to eqn (l.lla) and again B can be expressed in terms of microscopic quantities. As expected the introduction of partially slippery boundary conditions together with various models of viscoelasticity does not alter the long time behaviour of Co(t)which is given by eqn (1.4) (with D = 0). In fig. 1.3 molecular dynamics (dots) are compared with various hydrodynamic models with arbitrary slip. One curve corresponds to the purely viscous model in which the viscosity is chosen to give the best fit to experiment.Another curve corresponds to the purely viscous model in which the viscosity is computed using the Enskog viscosity. The last curve corresponds to the viscoelastic model of Theo-dosoupolu and Dahler in which all the parameters are computed by kinetic theory. In all of these models the coefficient of slip j3 is computed from the initial slope and the values are listed in the figure. It should be noted that in all cases the value of p is much closer to slip than stick. It should be noted that there are three distinct time regimes. First Co(t)decays exponentially with a rate determined from binary collisions. This is followed by a slower exponential decay with a renormalized rate constant and finally the Cw(t) decays with the asymptotic tail.O.OOI~-1 I I i 1 i I I I 0 4 8 12 16 20 24 30 time m mean collision times FIG.1.3 -A comparisonof the angular velocity correlation function for varioushydrodynamic models with molecular dynamics (dots). The reduced density is fi = 0.333. The slip parameter is denoted 8. The reduced viscosity is denoted $(qE = 0.855). Stick (best fit) $ = 0.103 B = 00 y = 0.0909. Slip (best fit) $ = 0.267 B = 1.97 y1 = y2. Slip (ab initiu) = 0.855 B = 0.465 y1 = 0.792 72 = 0.495. It is useful to determine when the best fit model reaches its asymptotic decay. This is given in fig. 1.4. BRUCE J. BERNE 53 Iin mean collision times 30 20’ 13 10 7 0.10 1 1 I 0.08 d 2 0.06 / Q3 0.04 0.02 0 1 0.005 0.010 i’2 in mean collision times FIG.1.4.-A plot of the best fit C&) against t -5/2.This curve looks fairly linear for times between 7 and 13 mean collision times but it is clear from the inset that this is not its time-asymptotic behaviour. The asymp- to tic behaviour becomes observable only beyond approximately 30 mean collision times. It is clear from this that although the preasymptotic behaviour may appear to behave like some power law it is dangerous to come to any conclusion based on data found for times shorter than 30 or 40 mean collision times. In concluding this section it is worth noting that for arbitrary slip boundary condi- tions and for all viscoelastic models the rotational diffusion coefficient is 2. COMPUTER EXPERIMENTS ON ROTATIONAL LONG TIME TAILS Given the remarkable agreement between the hydrodynamic prediction of the linear velocity correlation function [as embodied in eqn (1.5)] and various computer experiments it is natural to assume that the hydrodynamic prediction of the angular velocity correlation function [as embodied in eqn (1.3)] will likewise be in agreement with experiment.ROTATIONAL RELAXATION IN LIQUIDS O’Dell and Berne l5 have performed molecular dynamics experiments on fluids composed of 512 rough discs subject to periodic boundary conditions. Typical trajectories consisted of anywhere between 2 x lo6 and 5 x lo6collisions. The independent binary collision approximation for rough discs yields the short time Enskog exponential results for the linear and angular velocity correlation func- tions respectively,15 Cv(z)= exp -[::3 Cm(z) = exp -[1Y KIT where IC = 41/ma2,0is the diameter of a disc and z is in units of the mean collision time hereinafter abbreviated to rnct.FIG.2.1.-Log plots of the normalized linear velocity correlation functions computed for rough discs. Notice the arrow at approximately 22 rnct. This is the general area where we expect to see the sound “ kinks ”. K specifies the dimensionless moment of inertia and 2 specifies the density in units of the closest packed density. Fig. 2.1 shows a log plot of the normalized linear velocity autocorrelation functions (VCF) for three different svstems. For short times (<4 rnct) there is excellent agree- BRUCE J. BERNE ment with eqn (2.1).At longer times all three curves show a positive deviation from the exponential. At about 20 rnct there is a characteristic “ sound kink” in the VCF’s. This kink arises from the arrival of a sound wave from adjacent periodic cells and consequently illustrates gross interference effects due to the boundary conditions. We expect to find long time tails between 4 or 5 rnct up to a maximum of 20 rnct. For longer times the results are obscured by the boundary conditions. In fig. 2.2 corresponding plots are given for the angular velocity correlation func- tions (AVCF). These functions start deviating from eqn (2.2) at relatively short times when compared with the VCF’s. Since the same kind of hydrodynamic argu- ments apply to the long time tails in both Cu(z)and Cw(z),we expect that these tails will become manifest in the same time regimes.I .O c u \I I I I \. I I1 I 8‘ 12 16 20 24 28 T (in mean collision times 1 FIG.2.2.-Log plots of the normalized angular velocity correlation functions computed for rough discs. In fig. 2.3 C,(z) is plotted against 7-l. If there are “long time tails ” as predicted by eqn (1.5) for d = 2 these plots should be linear functions of z-l as z-l -f 0 as indeed they are. The coefficient a is determined from the slope and is presented in table 2.1. In fig. 2.4 Cm(z)is plotted against r2.If Cm(z)decays as r2at long times as predicted by eqn (1.3) with d = 2 these curves should also possess a definite linear region as they do. The coefficient aa as determined from the slope is presented in table 2.1.ROTATIONAL RELAXATION IN LIQUIDS Before discussing these results it is very important to note that as z -+ 00 the linear region Cm(z)does not extrapolate to zero. We think that this arises from the periodic boundary conditions. Our analysis of this effect is presented in a forthcoming manuscript. Had Alder and Wainwright not made their (N-1) correction their asymptotic linear velocity correlation functions would also extrapolate to a negative value. In their studies on smooth discs Dorfman and Cohen16have derived a formula TABLE 2.1.rALCULATED AND MEASURED VALUES OF THE DECAY COEFFICIENTS FOR THE LINEAR AND ANGULAR VELOCITY CORRELATION FUNCTIONS iia K a,(meas) a,(calc) a,(meas) a,(calc) runb 0.500 0.01 0.18 0.21 0.22 4.1 x 10-4 2.5 0.500 1.00 0.06 0.11 3.9 1.2 x lo-* 5.0 0.289 0.40 0.07 0.10 0.71 2.1 x 10-3 5.0 The density is measured relative to closest packed density.The number of collisions in the run in millions. 0.040L 16-5rcT 0.030-c. e c 0 .-Y 0 t 7 r 0 .-r K =0.4 -. 0.020-z aJ L 0 x c. 0 L Q) > L g 0.010-C .-d K = 1.0 n" = 0.500 n I I "0 0.1 0.2 I/T (in inverse mean collision times) FIG. 2.3.Plots of linear velocity correlation function against reciprocal time to display the long time tails on these functions for the fluid of rough discs. BRUCE J. BERNE from kinetic theory for the coefficient a,,which is identical to eqn (1.5) with the one crucial difference that v and D are replaced by their Enskog values vE and DB-a result that agrees well with the computer experiments of Wood and Erpenbecklg on hard discs.Unfortunately vE and DEare unknown for rough discs. O'Dell and Berne15 have performed this calculation. The values of a predicted on the basis of -0.020 16-5rcT ICT x h w u 2 C .-U V c 3 rc C .-L z 0.010-Q L V x e .-u 0 d Q) L b d 7 CI) 5 1/r2 ( in inverse mean collision times squared 1 FIG.2.4.-Plots of the angular velocity correlation function against reciprocal time squared to show the long time tails for rough disc fluids. this calculation are listed in column 4 of table 2.1.The agreement between GC determined from computer experiment and kinetic theory should be noted. The value of a determined from the kinetic theory is listed in column 5 of table 2.1. The most dramatic result of this study is that the predicted and "measured " values of ola differ by more than two orders of magnitude. The measured angular velocity long time tail is much larger than expected on the basis of all the theories. Were it not for the hydrodynamic theory presented in section 1 this result would be very disconcerting. Fig. 1.4 clearly shows that we are seeing a preasymptotic decay. Inthis connection Nady and Bernel' have completed a study of two dimen- sional fluids of ellipses containing 1296ellipses in which each pair of ellipses interacts according to a Lennard-Jones 12-6 potential with parameters e and Q which are dependent on the relative orientations of the ellipses.The form of E and 0 are given analytically by Berne and Pechukas.18 The dynamics are studied at different reduced ROTATIONAL RELAXATION IN LIQUIDS temperatures T* and different axial ratios a(=length/width). Cv(t)is shown in fig. 2.5. Again we see the sound kinks. The corresponding angular velocity correlation functions are presented in fig. 2.6. In addition the correlation function of the molecular orientation I 1 1 I L 010 0.20 0.30 I /t (inverse time in reduced units) FIG.2.5.-The linear velocity correlation function for Lennard-Jones ellipses for various reduced temperatures T*and various axial ratios a.(a) T*= 0.805 a = 2.0; (6) T*= 2.415 a = 1.0; (c) T*= 2.415 a = 1.3. is determined. This is shown for T* = 0.805 and a = 2.0in fig. 2.7. In these studies it is assumed that at long time C,(t) -+ aVt-l Cw(t)-t awt-2 and C,(t)3 CC~~-~. If the measured value of a is used to compute or it is once again found that the theoretical value of ucodiffers from the experimental value by two orders of magnitude-a result totally consistent with the foregoing. Thus we are observing a preasymptotic region here. In addition we see that the orientational correlation function behaves pre- asymptotically as does Co(t). In all of the foregoing studies it would appear that we are observing a preasymp- totic behaviour of the rotational correlation functions.This behaviour seems to last for an astoundingly long time and should have some important consequences. BRUCE J. BERNE I- t 1/t2 (inverse squaretime in reduced units) FIG.2.6.-The angular velocity correlation function for Lennard Jones ellipses for various reduced temperatures T*,and various axial ratios. (a) T* = 2.415 a = 1.3; (b) T*= 0.805 a = 2.0. Recently Desai and co-workers20 have calculated Cw(t)and C,(t)for the rough sphere fluid by mode-mode coupling theory and obtain reasonable agreement with our rough sphere calculations. In addition to this our own calculations indicate that periodic boundary conditions cannot be responsible for the large preasympto tic correlations. In closing we mention that in one run of 40 000 steps on a system of ellipses of low axial ratio (a = 1.3) and very high reduced temperature (T* = 2.415) there is clearly a region that looks asymptotic (see fig.2.8). The deviation from the exponen- tial is plotted against l/t2 in fig. 2.9. We have not yet completed data analysis so that fig. 2.9 should be regarded as tentative. Nevertheless it should be noted that the slope at long time is within a factor of 10 of what is predicted. ROTATIONAL RELAXATION IN LIQUIDS 0.040 I -0.030 - Y h-0.020-0.010 -I I ' O ' I 0.05 ' ' ' 0.10 ' ' ' 1 /t2 (inverse squaretime in reduced units) Fro.2.7.-The orientational correlation function defined in eqn (2.3) for Lennard-Jonesellipses at T*= 0.805 and a = 2.0. BRUCE J.BERNE S O .-c V t Y .e I= .-Y -L L 0 ....... . -0. .-0.OOll I 1 I I I I I f 1 I I I I I I 0 2.0 4.0 6.0 8.0 10.0 12.0 14.0 reduced time It/t,J FIG.2.8.-CV(t) and Cw(t)for a fluid of ellipses of axial ratio a = 1.3 and reduced temperature T*= 2.415. At = 5 X 10-3t,. to = (ma"/&)+. *.* ROTATIONAL RELAXATION IN LIQUIDS t units of t/f c0.008 - 0.006- a .. 1 - *. 0.004-a0 -0.002--0 lll’1llllllll roughly indicates small amplitude of long time tail in Co(t). T* = 2.415 a = 1.3 p/po-= 0.3 N = 1296. We have benefited from several very useful conversations with Prof. R. Kapral and C. Hynes. These investigators have independently developed a hydrodynamic model.B. J. Alder D. M. Gass and T. E. Wainwright J. Chem. Phys. 1970 53 3813. C. M. Hu and R. Zwanzig J. Chem. Phys. 1974,60,4354. G. K. Youagren and A. Acrivos J. Chem. Phys. 1975 63 3846. D. R. Bauer J. I. Brauman and R. Pecora J. Amer. Chem. Sac. 1974,96,6840. R. Zwanzig and M. Bixon Phys. Rev. A 1970,2,2005. B. J. Alder and T. Wainwright Phys. Rev. Letters 1967 18 968. L. Verlet D. Levesque and G. Weiss personal communication. N. K. Ailawadi and B. J. Berne IUPAP Conference on Statistical Physics Chicago March 1971. B. J. Berne J. Chem. Phys. 1972,56,2164. lo B. J. Berne unpublished manuscript. l1 F. Garisto and R. Kapral Phys. Rev. A 1974,10,309; Y. Pomeau and J. Weber Phys. Rev. A 1973 8 1422; T. Keyes and I. Oppenheim Physica 1974,75 583.l2 B. Hills and J. M. Deutch Physica 1976 in press. l3 J. Montgomery and B. J. Berne J. Chem. Phys. 1976 (2 papers in press). l4 J. O’Deil and B. J. Berne J. Chem. Phys. 1975 63 2376. l5 J. O’Dell and B. J. Berne The Persistence of Linear and Angular Velocity in a Fluid of Rough Discs Manuscript in preparation Sept. 1976. See also J. O’Dell Ph.D. Dissertation (Columbia University 1976). l6 J. R. Dorfman and E. G. D. Cohen Phys. Rev. Letters 1970 25 1257. l7 L. Nady and B. J. Berne Persistence Effects in Anisotropic Fluids manuscript in preparation Sept. 1976. l8 B. J. Berne and P. Pechukas J. Chem. Phys. 1972,56,4213. l9 See for example W. W. Wood Fundamental Problems in Statistical Mechanics IIZ,ed. E. G. D. Cohen (North Holland Amsterdam 1975).2o J. R.Mehaffrey R. Desai and R Kapral J. Chem. Phys. 1977 66 1665.