Volume 163, number 4,5
CHEMICAL PHYSICS LETTERS
NON-MAXWELL VELOCITY DISTRIBUTIONS
I5 November 1989
IN EQUILIBRATED FLUIDS
S.-B. ZHU, J. LEE and G.W. ROBINSON Picosecond and Quantum Radiation Laboratory. Texas Tech University, P.O. Box 4260. Lubbock, TX 79409, USA Received I8 August 1989
In this paper, using supercomputer molecular dynamics methods, we further demonstrate that the velocity distribution of a low-mass classical particle subjected to a strong force field in a heat bath of heavy particles is non-Maxwellian. When the system reaches equilibrium, the non-Maxwellian distribution persists for velocity components along the field direction.
1. Introduction The Maxwellian velocity distribution was initially derived for spatially homogeneous gases [ I]. The Maxwell distribution remains valid for a thermally isolated ensemble where the position and momentum of each particle form a pair of independent stochastic variables. This distribution applies even in the case where there may exist spatial inhomogeneities caused by an external conservative force field. In a condensed phase, the interactions are generally of the many-body type. .4ccording to the Langevin theory, the total force on a particle can be divided [ 21 into a random part and a systematic part, the latter giving rise to the friction term. If the process is non-Markovian, the friction can be expressed as an integration of the memory kernel. Because of the “heat bath” friction, the system particle experiences a modified potential that is different from its gas phase value [ 3 1. The corresponding force is generally not conservative. As a matter of fact, the modified quantity may no longer be a genuine potential, since it depends not only on the position but also on the history of the particle. For nonlinear interactions, the memory kernel depends in a complicated way [4] on both the position and the momentum of the particle on interest (the “system particle”). We fell that it is the presence of these nonlinear effects that destroys ordinary intuition about this problem. Because of these types of complications, the effective Hamiltonian of the system particle cannot in general simply be written as a sum of a potential en328
ergy, which is a function of position only, and a kinetic energy, which is a function of momentum only [ 51. In other words, the ensemble may no longer be a canonical one. It has never been clear whether such systems possess the standard Maxwellian velocity distribution, though most authors have assumed that they do. Understanding this phenomenon is certainly of significance in many dynamical problems, including methods used for velocity sampling in approximate MD [6] and Monte Carlo calculations [ 71, or in the kinetics study of H- or D-atom barrier crossing problems in solids or liquids #I. Recent full-scale computer molecular dynamics (MD) studies [ 10, 1 1 ] of Kramer+type [ 12,131 potential-barrier crossings in condensed phase cnvironments had suggested earlier that local noncanonical velocity distributions might sometimes persist for certain types of ultrafast dynamical processes. The question posed here is: “Just how fundamental and widespread are these deviations?” In other words: “Do they occur even for particles that do not possess the added complications of internal vibrations and structure?’ Because of the current interest in the dynamics of low-mass atoms interacting strongly with a heavy particle environment [ 8,9], MD computations were recently performed on an ensemble of deuterium-like atoms interacting strongly with a twoxl For example, the phenomena reported here may play a role in H-atom diffusion processes in metals or semiconductor materials [ 81and possibly in the recently reported observation of cold nuclear fusion [ 91.
0 009-2614/89/$ 03.50 0 Elsevier Science Publishers B.V. (North-Holland Physics Publishing Division )
Volume 163,number 4,5
CHEMICALPHYSICSLETTERS
dimensional harmonically bound palladium-like lattice [ 14,151. Again, this system exhibited a nonMaxwellian distribution for the low-mass particles. In the present paper, our interest is in determining whether or not this unexpected phenomenon occurs also in a fluid. A further motive is to convince ourselves and others of the reality of the phenomenon by employing a totally different system than in any of our past work. It would be inconceivable that so many different systems, having different starting conditions, etc., could show the same behavior if it were not real. The only similarity between all these systems is the presence of low-mass particles (or a particle) interacting with a heat bath of heavy particles: the low-mass particle(s) also being subjected to very strong force fields.
2. Methods System 1. We begin this investigation by studying a pure Lennard-Jones (LJ) fluid with pair-wise additive interparticle potentials u(r) =4~ [ (a/r) ‘* (c~/r)~]. Each particle is also subjected to an external potential applied along the x direction. This sample is composed of 108 identical particles (system I), and is extended in the three Cartesian directions through periodic boundary conditions. Standard NVE MD methods [ 161 are used with an integration time step At=O.O05om, where M is the mass of the LJ particle. Besides the LJ interactions, each of these particles experiences an additional time-invariable force caused by the applied potential. This potential is assumed to have the following form,
U(x) = Liucos(4Rx/l)
)
(1)
where L is the length of the cubic sample, and U0 is taken to be 10 k,T. The period 4L in system I is about 2.56 times the LJ diameter u. Since U(x) oscillates across the sample in a spatial sense only, there is no tendency for the particles to drift from one side of the box to the other. The particles thus spend most of the time in the vicinity of the two potential minima near x= !L and x= fL, and occasionally cross the barriers near x=0, &, L. For such an ensemble, the external force basically gives rise to a dissymmetry in the Cartesian directions and to a cooperative many-body response of
17 November 1989
the particles. In spite of these features, it was found that velocity autocorrelation functions for components parallel to and orthogonal to the force direction differ only slightly, even though a (free-atom) potential barrier of 20k,Tis present. Of course, just as for the two-dimensional harmonically bound heavy-atom lattice [ 14,151, the effective condensed phase potential barrier is much lower than the freeatom barrier height. In addition, the minima are shifted by = & &L because of the strong particleparticle repulsions from the LJ forces. One of the disadvantages of computer simulation is the limitation on the number of particles that can be used in Newton’s equations. To avoid this problem, a long time average for a steady-state replaces the ensemble average [ 16 I. The rapid evolution from a nonequilibrium to an equilibrium state has been well known in MD calculations for many years #2. Using a 1.2 x lOSAt run after reaching equilibrium, we obtain perfect Maxwellian velocity distributions for system I at any value of x, including the region near the barrier top. These Maxwell distributions v&e found to set in very early ( < 104At) and are independent of whether they are computed along or normal to the external field direction. The probability of any of the 108 particles being near the barrier top, 06 Ix+!&,
329
Volume 163,number 43
CHEMICALPHYSICSLETTERS
depends on the specific position of the system particle along the force directions. For example, when the system particle is near a minimum of the applied potential, high frequency oscillations appear in the velocity autocorrelation function of components along the force direction, but not for the two transverse components. The above situation creates a circumstance similar to the one described by the generalized Langevin equation. Three forces affect the motion of the system particle. The first one, the random force from neighbouring bath particles, provides the main actuation for disorder which, alone, would eventually lead to a canonical distribution. The other two forces, the dissipative bath friction and the conservative force from the external potential, belong to systematic variables that tend to increase the order. In nonlinear systems, fluctuations can also produce order [ 181. Short-time information can survive to affect long-time behavior, even in the presence of intense stochastic forces [ 191. The resulting distribution is determined through the competition of these three effects. If the tendency to produce order is sufiiciently strong, a departure from the canonical distribution may become possible. To increase the effects of the barrier force, and thus to magnify the tendency towards order so that a measurable effect is created, we reduce the period of the external potential in eq. ( 1) to 0.02L. In order to improve the statistic, we take U. in system II to be 2k& giving a 4kBT free-atom energy from minimum to maximum. This is equivalent to a barrier height of 6.06~,,+,, where emM=fiW denotes the LJ parameter of interaction of the system particle with a heat bath particle. Various (4.92-7.32) x 106At runs indicate nearly canonical Maxwellian distributions at any positions, including the vicinity of the barrier maxima, O< Ix+O.O2nLI <0.0025L(n is an integer), for motion of the low-mass particles transverse to the external field. See fig. la. Just as for system I, these Maxwellian distributions set in very rapidly compared with the time scale of the entire run. However, as shown in fig. 1b, there is a severe distortion from the Maxwellian distribution for velocity components of the system particle along the barrier force direction. This distribution also sets in very quickly and remains intact for the remainder of the run. 330
17 November 1989
J.O Fig. 1,(a) Distribution functions (normalized to I ) for velocity componentsof the systemparticletransverseto the externalforce. The reduced velocityis defined as v/ (v*) I’*,where v2= us + u,‘. The curve denotes the Maxwelliandistribution: C,vexp( -v2/ (v*) ) dv, the symbolsrepresent the simulation results,with circles denotingthe total ensembleaverageddistributionand crosses the distribution near the barrier top. Data averaged over 4.92x 10%. (b) Distribution functions for the longitudinalvelocity.The reduced velocityis defined as u,/( v: ) i/z. The curve denotestheMa.xwelliandistribution:C, exp( - vf/(v:)) dvi the crosses representdata averagedover 3.72~ 106Af,the circlesover 7.32x 1064f.
The distribution in fig. lb is seen to be fairly flat for v, 5 (0:) ‘j2, where (v: ) denotes the mean square velocity of the system particle. It then falls off when v,/ ( v: ) ‘I2 increases further, but decays more slowly than Maxwellian for high velocities. This latter feature, the “high tail” of the distribution, was also found in our earlier work [ 14,15 1. This behavior could be crucial lo the understanding of kinetic phenomena for low-mass atoms in a heavy atom environment [ 8,9], since there may be many times the number of hot low-mass particles than given by the Maxwellian distribution. The results indicate a significant noncanonical distribution even in the barrier minimum region. However, the largest deviations appear on the way
Volume 163,number 4,5
CHEMICALPHYSICSLETTERS
to climbing to the barrier top, where the barrier force is the strongest and the low-mass particle acceleration is the greatest. Kinetic energy contour maps can most easily be visualized in two-dimensional problems. These are depicted in pending publications [ 14,151. In the absence of the external barrier force, the distribution was found to recover its Maxwellian form in accordance with earlier findings [ 10,l I 1.The choice of t, does not affect the results sensitively. For example, an increase (or decrease) of E, by a factor of 16 ( em,uby a factor of four) leads to similar distributions. Of course, as cm-O, we approach the free oscillator velocity distribution which has a high peak near the cut-off maximum velocity, then suddenly drops towards zero beyond that point 113.
3. Effect of nonlinearities In 1962, Lebowitz and Rubin [ 211 studies the behavior of a light particle in a fluid from a dynamical point of view. When a constant unidirectional external force was allowed to act on the system particle, these authors obtained a stationary noncanonical distribution representing a balance between the effect of the acceleration caused by the external force and scattering by the bath particles. The deviation from the canonical distribution depended on the strength of the force and the mass ratio of the system particles to a bath particles. In the present problem, the force is not constant across the sample. The external force accelerates the system particle at a certain position, but decelerates it at some other position. Any deviations from normal velocity behavior would tend to cancelif the sys-
tem is linear, and one would obtain a canonical distribution. However, in the presence of nonlinearities, i.e. when the perturbation does not linearly depend on the force, or when nonlinear state-dependent friction exists, this cancellation is not possible. The history of the system particle comes into play, giving different velocity behavior, for example, depending upon whether the system particle is moving towards the barrier top or away from it. When these nonlinear elements are large relative to the stochastic ele#’Similar to the classical spatial distribution of a free one-dimensional harmonic oscillator.See, for example,ref. [ 201.
17 November 1989
ments, large deviations from the Maxwellian distribution seem to be possible. In fact, additional computations [ 221 indicate that distributions spanning the entire range from pure Maxwellian to free oscillator [ 201 can be obtained.
4. Conclusions Is it still a possibility that the noncanonical distributions that we have now observed in so many different systems are brought about by an unphysical reason such as inadequate run times? Our answer is no. The data obtained for system II have been averaged over 7.32~ 106Ait,which, as far as we are aware, is considerably longer than in any previously published work from any other laboratory. It would be expected that after reaching equilibrium all the distributions found here, Maxwellian or not, would set in well within 106At. Notice, for example, in fig. la that a near Maxwellian velocity distribution for the system II low-mass particle in the vicinity of the pbtential maxima is observed along the y and z directions, though taken from only 4.92X lo6 x 1.47% = 7.23 x 1O4samples. (The percentage here is the measured probability for a particle to lie within 15% of the barrier maximum.) The non-Maxwellian distribution in fig. lb has been observed to remain virtually constant from 1 X 106At through the entire (3172-7.32) x 106At run. Finally, the probability of the system particle having positive velocity components along the field direction is the same (within 0.05%) as those having negative velocity components. This is the expected result from the principle of detailed balance for the potential used. It is difficult to see how all these findings in this and other diverse systems could have occurred if the run times or statistics were inadequate. The results seem also to tell us that the noncanonieal distribution is essentially not caused by rotations or translations of the “reaction coordinate” in the space-fixed frame of ref. [ lo]. It is more likely a consequence of nonlinear and non-Markovian processes combined with order-disorder competition. When the low-mass system particle moves rapidly under the influence of the strong systematic force, it can never reach perfect equilibrium with neighboring slow moving bath particles. In other words, the 331
Volume163,number4,s
CHEMICAL PHYSKCS LETTERS
interface between the system and its heat bath is not one of total disorder. The position and the momentum of the system particle lose their independence, and are no longer a pair of canonical variables. In conclusion, by reference to our earlier work [ 10. I 1,14,15 1, and by comparing the results for systems I and II here, we find that, ifa system particle has a low or comparable mass relative to the buth particles, and jiirther ifthis low-mass particle is subjected to a strong force, the particle may then show a severe departure porn a canonical distribution. We believe this may be an important characteristic of many materials in condensed phases.
Acknowledgement Financial support at the PQRL has been shared jointly by the Robert A. Welch Foundation (D-0005, 50% and D-1094, S%), the National Science Foundation (CHE8611381, 39%) and the State of Texas Advanced Research Program (1306, 6%). Computer time was furnished by the Texas Tech University Academic Computing Center and the Pittsburgh Supercomputing Center. Helpful discussions with Dr. S. Singh and Professor R.K. Pathria are gratefully acknowledged.
References [ 1] F.W. Sears, An introduction to thermodynamics, the kinetic theory of gases and statistical mechanics, 2nd Ed. (AddisonWesley, Reading, 1953) ch. 12. [2] H. Mori, Progr. Theoret. Phys. 28 (1962) 763.
332
17 November 1989
[3] S.-B. Zhu and G.W. Robinson, Entropy effects in liquid phase isomerizations, Proceedings of the 3rd International Supercomputing Institute., Vol. 1 ( 1988) 300. [4] S.-B. Zhu, J. Lee and G.W. Robinson, .I. Chem. Phys. 88 (1988) 7088. [ 5] S.-B. Zhu, S. Singh and G.W. Robinson, Phys. Rev. A 40 (1989) 1109. [6] T. Yamamoto, J. Chem. Phys. 33 (1960)28 1. [7] N. Metropolis, A.W. Metropolis, M.N. Rosenbluth, A.H. Teller and E. Teller, J. Chem. Phys. 21 (1953) 1087. [S] S. Estreicher, Phys. Rev. A 36 (1987) 9122; S.J. Pearton, J.W. Corbett and T.S. Shi, Appl. Phys. A 43 (1987) 153. [ 4] S.E. Jones, E.P. Palmer, J.B. Czlrr, D.L. Decker, G.L. Jensen, J.M. Thome, SF. Taylor and J. Rafelski, Nature 338 ( 1989) 737; M. Fleischmann, S. Pons and M. Hawkins, J. Electroanal. Chem. 261 (19B9) 301. [IO] S.-B. Zhu and G.W. Robinson, Chem. Phys. Letters 153 (1988) 539. (I I ] S.-B. Zhu and G.W. Robinson, J. Phys. Chem. 93 (1989) 164. [ 121 HA Kramers, Physica 7 ( 1940) 284. [ 131 S. Chandrasekhar, Rev. Mod. Phys. 15 ( 1943) 1. [ 141 S.-B. Zhu, J. Lee and G.W. Robinson, Chem. Phys. Letters 161 (1989) 249. [IS] S.-B. Zhu, J. Lee and G.W. Robinson, Non-Maxwell velocity distributions in inhomogeneous materials, Proceeding Workshop on Cold Fusion Phenomena, Sante Fe, NM (May 23-25 1989) (Plenum Press, New York), submitted for publication. [ 161 L. Verlet, Phys. Rev. 159 ( 1967) 98. 1171 P.L. Fehder, J. Chem. Phys. 50 (1969) 2617. [ 181M. Suzuki, Advan. Chem. Phys. 46 (1981) 195. [ 191 W.T. Coffey, M.W. Evans and P. Grigolini, Molecular diffusion and spectra (Wiley-lnterscience, New York, 1984) chs. 7 and 8. [201 A. Yariv, Quantum electronics, 2nd Ed. (Wiley, New York, 1975) pp. 24-25. 1211 J.L. Lebowitzand E. Rubin, Phys. Rev. 131 (1963) 2381. 1221 S.-B. Zhu, J. Lee and G.W. Robinson, unpublished work.