- Email: [email protected]
Long-time tails in two-dimensional fluids by molecular dynamics
PHYSICA ELSEVIER
Physica A 240 (1997) 268-276
Long-time tails in two-dimensional fluids by molecular dynamics Mauro Ferrario a,*, Antonino Fiorino b, Giovanni Ciccotti c alstituto Nazionale per la Fisica della Materia and Dipartimento di Fisica, Universitdt di Modena, Italy b Istituto Nazionale per la Fisica della Materia and Dipartimento di Fisica, Universith di Messina, Italy c Istituto Nazionale per la Fisica della Materia and Dipartimento di Fisica, Universitgt di Roma "La Sapienza", Italy
Abstract We report on molecular dynamics simulation of long-time tails in the velocity and stress autocorrelation functions of a dense two-dimensional fluid. Large systems of the order of hundred thousand particles have been investigated, performing canonical averages over an ensemble of trajectories generated on a parallel computer. We find the well-known t-~ decay for the velocity autocorrelation function at two different densities of the fluid, together with a faster than linear time dependence for the mean-square displacement at long times. Although there are indications of an asymptotically faster decay, the data are not precise enough to discriminate whether the decay is in agreement with the (t Ix/i-n--/)-J prediction of consistent mode-coupling theory or it is due to finite size effects. No evidence, within the statistical errors, is found for a long-time tail in the stress autocorrelation function. This finding is in agreement with recent NEMD results [Hoover et al., Phys. Rev. E 51 (1995) 273; Gravina et al., Phys. Rev. E 52 (1995) 6123], who find an analytical dependence of the shear viscosity upon the shear rate with no evidence for divergence in the Green-Kubo value.
1. Introduction The validity of the Green-Kubo formulas relating transport coefficients to equilibrium fluctuations has been a subject arousing considerable controversy. After the simulation o f hard disks by Alder and Wainwright [1], who found the presence o f the long-time tail in the velocity autoeorrelation function beyond the expected exponential decay on the scale o f the mean collision time, despite many efforts [ 2 - 6 ] , no qualitatively new results have been obtained to date by molecular dynamics. Only recently, by lattice-gas cellular automata [ 7 - 9 ] , Frenkel and coworkers have been able to identify the detailed behavior o f long-time tails. * Corresponding author. Tel.: + 39.59.586059; Fax: +39.59.367 488; e-mail: [email protected]. 0378-4371/97/$17.00 Copyright (~) 1997 Elsevier Science B.V. All rights reserved PII S 0 3 7 8 - 4 3 7 1 ( 9 7 ) 0 0 1 50-7
M. Ferrario et al./Physica A 240 (1997) 268-276
269
On the theoretical side an approach, based on mode-coupling theory, predicts [10], in two dimensions, a power-law decay of the form (v(O)v(t)) o( [(D + v)t] -1 ,
(1)
where D and v are assumed to be constant, and are, respectively, the self-diffusion coefficient and the kinematic viscosity, related to the shear viscosity q by ~/= mpv. The above result, not consistent in two dimensions as it leads to the In t divergence of the, assumed constant, coefficients D and v, has been improved by carrying out a selfconsistent calculation. The self-consistent approach, using time-dependent coefficients [11,12] D(t) and v(t), yields the asymptotic behavior [13]
(v(O)v(t)) cx [(D(t) + v(t))t] -i ~ (t Ix/~-~n t) - I ,
(2)
which gives D(t) cx lx/~7. A similar behavior holds for the stress autocorrelation function, and leads to a slow divergence in two-dimension also for the viscosity coefficient 11. The aim of the present study is to calculate the velocity and the stress autocorrelation functions of a fluid system with a number of particles sufficiently large for the asymptotic regime to be reached and with as much statistics as possible to make an accurate test of these predictions even at high density.
2. Details of the computation We have chosen to study a 2D s o . d i s k fluid, whose particles interact by a pair potential of the form
u(r) = 4~(a/r) 12
(3)
at two state points characterized by the same temperature, kBT/e = 1 and a mass density of pa2= 0.69 and po'2= 0.96, respectively, the latter close to the solid transition [14]. Extensive MD-simulations were performed with a cut-off radius r~. = 1.5 and a time step 6t = 0.005, both values are in reduced m, a, e units, with a corresponding time-unit to = a(m/e) 1/2. The most accurate calculations were performed for a system of 65 536 (64 K) particles and in order to explore system size dependence also systems of 36864 (36K) and 16384 (16K) were investigated. For smaller systems boundary condition effects, related to the recurrence time, interfere severely with the slow, long-time, power-law decay. According to the value of the velocity of sound c = 9.08 calculated by Gravina et al. [14], for the high-density state point the recurrence time is zL =L/c = 28.7 for the 6 4 K particle system, zL = 21.6 for the 3 6 K particle system and zL = 14.4 for the 16 K particle system. At the lower density we expect these times to be even longer, for the combined effects of the increase of the box length and the decrease of the density. The equations of motion were integrated using the velocity-Verier algorithm with on the fly calculation of the autocorrelation functions of the particle velocity Vi(t) and
270
M. Ferrario et al. I Physica A 240 (1997) 2 6 8 ~ 7 6
of the shear stress ~xy = _ _
....
(4)
i:l
j~i
At the same time the mean-square displacement was also directly computed
Ar2(t) =
1
N
~ Z
<[ri(t) -- ri(O)]2> .
(5)
i=1
A large statistical sample was produced using a massively parallel computer, a Cray T3D with 128 processor, with each node computing an independent system trajectory. Starting from a unique initial condition each system was carefully equilibrated by first running it at higher temperatures, with a different, randomly chosen, value of the temperature for each node and subsequently, each system was thermalized applying a sequence of velocity randomizations, using an independent random number generator on each node. During equilibration, trajectories were evoluted according to Nos6-Hoover constanttemperature equations of motion [15]: /'i = P__L,
P~
Pi = Fi -- --~'P',
N
P~ // = -Q-'
t5, = Z
mi
P~ -- -
2NksT,
(6)
i= 1 mi
using a modified velocity-Verlet reversible scheme of integration derived using the Trotter-Suzuki factorization of the Liouville operator [16]. Using a five-term decomposition of the Liouville operator iL, with Ui(z) = exp{iLiz), the approximate discrete time propagator G(t) G(T)-.~ U5(2)U4(2)U3(2)U2(2)Ul(T)U2(2)U3(2)U4(2)U
5
(2)
(7)
is constructed with the following break-up
iL = iLl + iL2 + iL3 + iL4 + iL5 =
F,.N+
.
.
.
.
+ p~. @.
(8)
The resulting algorithm has the following structure, for z = At:
ri(At)
= ri(O) + At ~i(At/2),
r~(At) = {~(At/2) + ~--tmiF~[{r(At)}])exp { - 2 ~(At/2)), q(At) = q(O) + Ate(At/2), At
O(At) = ?7(At~2) ÷ 5-QF~[(~(3t))],
(9)
M. Ferrarioet al./PhysicaA 240 (1997) 268-276 where one first evaluates the half-step values
~(At/2) :
271
¢I(At/2) and ~i(At/2):
r~(0) + ~A--~-t F,[{/'(0)}],
~i(At/2):fi(O)exp {-A---~~(At/2)} + 2A--~tmiFi[{r(O)}]
(10)
and the forces in the above equations are given by Fi[{r(t)}] -
c~U[{r}] , ~ri {r(t)}
F,[{f(t)}] = p~[{/'(t)}] .
(11)
Standard microcanical equations of motion were integrated thereafter and the choice of the velocity-Verlet scheme (or of the modified algorithm used for Nos~-Hoover dynamics) insured the minimum memory usage, since only one 2D position and one 2D velocity per particle needs to be kept in memory at any time. This allowed us to fit large systems on the limited RAM of each node. In fact, the size of the system was primarily limited by memory considerations related to the calculation of the autocorrelation functions along the run, rather than by the computation of the dynamical trajectories. Although systems of more than 100 thousand particles were tested, most calculations were for a system with 16K, 36K and 64K, where K stands for 1024 multiples, this choice being dictated by architectural I/O considerations. Reported results for the 65 536 particle system correspond to a global statistical average over almost 18 millions steps for the low density and more than 30 millions for the high-density point.
3. Results The autocorrelation functions have been calculated for the 64K particles system up to a time t = 40.96 for both state points. At the lower density the velocity autocorrelation function is positive definite in all the time region considered here. It exhibits a shorttime exponential decay, with also a relatively small minimum at t = 0.5 followed by the long-time t -1 regime which approximately starts at time t = 1 and remains very evident all along up to t ~ 4 1 as it is shown in Fig. 1. At the higher density, see Fig. 2, the velocity autocorrelation function has a much more oscillatory behavior at short times and a number of negative-valued minima are present. At around time t = 1 the function becomes positive again and start decaying slowly with a clear t - I power law. For this density the long-time tail is one order of magnitude smaller than for the previous case. Also around time t = 30, that is about the estimated recurrence time, there is a clear evidence of possible interference effects due to periodic boundary conditions. A more careful analysis of the results for the lower density shows some evidence for a somewhat faster than t - t behavior for the larger times. In fact, a least-squares fit of the long-time tail, shown in Fig. 3 gives an exponent o f - 1 within 10 -3, more
272
M. Ferrario et al./ Physica A 240 (1997) 268-276
ill 0.5
0
0
1
2
time
3
4
Fig. 1. Velocity autocorrelation function, p = 0.69. In the inset the long-time behavior is shown using a logarithm scale for both axis. 2
0
0.2
0.4 time 0.6
0.8
Fig. 2. Velocity autocorrelation function, p = 0.96. In the inset the long-time behavior is shown using a logarithm scale for both axis. precisely the exponent is 2 = - 1 . 0 0 0 2 , when the time region between t = 2 and t = 12 is taken into account, while a significantly larger deviation, 2 = - 1.05, is measured when using the far region, from t = 20 up to t = 40.96 (for the entire region a 2% deviation is obtained). Even with the high accuracy and extended range of the present data we feel that it is not possible to discriminate between a genuine faster than t-J behavior and finite-size effects. However, when in Fig. 4 we make a comparison of the velocity autocorrelation functions obtained with different system size for the low-density system, the results practically coincide in the overlapping time region. Notice that the results for the two
273
M. Ferrario et al./Physica A 240 (1997) 268-276
0.004 0.03 0.003
0.01 0.008
0.002
0.006 I
I
2
3
I
I
I
I
I
I ]
20
4 5 678910
30
40
Fig. 3. Power-law fit to the velocity autocorrelation function. The two used time region are explicitly shown, 1.2 0.1
0.8 O.Ol
~
6
K
0.6 0.4
0.001 1
10
0.2
0
l
2 time 3
4
5
Fig. 4. System size dependence for the normalized velocity autocorrelation function, p = 0.69. Curves have been appropriately shifted for a clearer presentation. smaller systems, 16K and 3 6 K particles, are averages performed on a single, long, microcanonical trajectory produced in several weeks on a dedicated workstation. The direct calculation of the mean-square displacements gives results which agree with the slow decay of the velocity autocorrelation function, that is the mean-square displacement grows faster than linearly with time, indeed proportional to In t, as it is clearly shown in Fig. 5 for the system at the lower density. The small deviation, expected according to the previous discussion, is not noticeable on this scale. The situation for the shear stress autocorrelation function is less conclusive. The collective nature of the observable and the overwhelming importance of the potential contribution makes it difficult to obtain high-precision data for the long-time tail, which is expected to be of kinetic origin, and it is not possible to discern a region
274
M. Ferrario et al. IPhysica A 240 (1997) 268-276
0.8 p=.69 0.6
0.4
0.2
0.01
0.1
time 1
10
Fig. 5. Ratio of the mean-square displacement over the time Ar2(t)/t.
8# 105
6 105
104
4 105
2111
103
,
0
,
,
,
I
1
,,,I
,
,
,~-~ 1
,
0.01 ,
0.1 ,
I
I
2
3
time
,
,
,
,
I
,
,
,
,
4
Fig. 6. Shear-stress autocorrelation function, p = 0.69. In the inset a logarithm scale is used for both axis, where a clear power-law decay emerges from the statistical noise. It has to be noticed, however, that, while in the velocity case the magnitude o f the long-time tail in the normalized autocorrelation function is some 3% in the low-density case and 0.2% at the higher density, for the stress autocorrelation functions we do not see any such feature at long times within a statistical noise o f 0.4% at lower density and 0.2% at higher density. The autocorrelation functions for the two systems are shown in Figs. 6 and 7. W h e n the shear-stress autocorrelation functions are integrated over time a welldefined plateau is obtained from which we can extract the G r e e n - K u b o values for the shear viscosity r / = 1.05 + 0.10 at p = 0 . 6 9 and r / = 5.91 4-0.04 at p = 0 . 9 6 , the latter in agreement with the result o f Gravina et al. [14].
M. Ferrario et al./Physica A 240 (1997) 268-276
275
2 106 106
lOS
~k~.
i io 6
0 0
1
2
time
3
4
Fig. 7. Shear-stress autocorrelation function, p = 0.96. In the inset a logarithm scale is used for both axis. 4. Conclusions
Our molecular dynamics calculations are consistent with a high degree of precision with the expected t -1 behavior for the long-time tail of the velocity-autocorrelation function of a two-dimensional fluid. Although some evidence is given for a somewhat faster decay the problem of the nature and the origin of this correction cannot be assessed with the present data. On the contrary, despite the effort, it has not been possible to establish the presence of the long-time tail for the stress-autocorrelation function. There are various possible reasons for that: (i) lack of statistics; (ii) the dominance of the potential contribution to the shear stress, orders of magnitude larger than the kinetic one; (iii) destructive interference of the vortices associated with the individual particle velocity on the overall time memory of the collective observable.
Acknowledgements
We wish to acknowledge enlightening discussions with Daan Frenkel and a critical reading of the manuscript by Ray Kapral. Computer time on the Cray-T3D at the National Supercomputing Center, CINECA, was made available through an initial grant allocation by Cineca and subsequently within the INFM parallel computing initiative.
References
[1] B.J. Alder, T.E. Wainwright, Phys. Rev. A 1 (1970) 18. [2] W.W. Wood, in: E.G.D. Cohen (Ed.), Proceedings of the Summer School on Fundamental Problems in Statistical Mechanics III, North-Holland, Amsterdam, 1975. [3] W.E. Alley, B.J. Alder, Phys. Rev. Lett. 43 (1979) 653. [4] D. Levesque, W.T. Ashurst, Phys. Rev. Lett. 33 (1974) 277.
276 [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16]
M. Ferrario et al./Physica A 240 (1997) 268-276
J.J. Erpenbeck, W.W. Wood, Phys. Rev. A 26 (1982) 1648. J.J. Erpenbeck, W.W. Wood, Phys. Rev. A 32 (1985) 412. M.A. van der Hoef, D. Frenkel, Phys. Rev. A 41 (1990) 4277. M.A. van der Hoef, D. Frenkel, Phys. Rev. Lett. 66 (1991) 1591. M.A. van der Hoef, D. Frenkel, A.J. Ladd, Phys. Rev. Lett. 67 (1991) 3459. For a review, see: Y. Pomeau, P. Rrsibois, Phys. Rep. 19 (1975) 63; see also: M.H. Ernst, B. Cichocki, J.R. Dorfman, J. Sharma, H. van Beijeren, J. Stat. Phys. 18 (1978) 237. B.J. Alder, T.E. Wainwright, D. Gass, Phys. Rev. A 4 (1971) 233. R. Zwanzig, Physica A 118 (1983) 427. M.A. van der Hoef, Simulation Studies of Diffusion in Lattice-Gas Fluids and Colloids, Ph.D. thesis, Utrecht, 1992. D. Gravina, G. Ciccotti, B.L. Holian, Phys. Rev. E 52 (1995) 6123. W.G. Hoover, Phys. Rev. A 31 (1985) 1695. M. Tuckermann, B.J. Berne, G.J. Martyna, J. Chem. Phys. 97 (1992) 1990.
Physica A 240 (1997) 268-276
Long-time tails in two-dimensional fluids by molecular dynamics Mauro Ferrario a,*, Antonino Fiorino b, Giovanni Ciccotti c alstituto Nazionale per la Fisica della Materia and Dipartimento di Fisica, Universitdt di Modena, Italy b Istituto Nazionale per la Fisica della Materia and Dipartimento di Fisica, Universith di Messina, Italy c Istituto Nazionale per la Fisica della Materia and Dipartimento di Fisica, Universitgt di Roma "La Sapienza", Italy
Abstract We report on molecular dynamics simulation of long-time tails in the velocity and stress autocorrelation functions of a dense two-dimensional fluid. Large systems of the order of hundred thousand particles have been investigated, performing canonical averages over an ensemble of trajectories generated on a parallel computer. We find the well-known t-~ decay for the velocity autocorrelation function at two different densities of the fluid, together with a faster than linear time dependence for the mean-square displacement at long times. Although there are indications of an asymptotically faster decay, the data are not precise enough to discriminate whether the decay is in agreement with the (t Ix/i-n--/)-J prediction of consistent mode-coupling theory or it is due to finite size effects. No evidence, within the statistical errors, is found for a long-time tail in the stress autocorrelation function. This finding is in agreement with recent NEMD results [Hoover et al., Phys. Rev. E 51 (1995) 273; Gravina et al., Phys. Rev. E 52 (1995) 6123], who find an analytical dependence of the shear viscosity upon the shear rate with no evidence for divergence in the Green-Kubo value.
1. Introduction The validity of the Green-Kubo formulas relating transport coefficients to equilibrium fluctuations has been a subject arousing considerable controversy. After the simulation o f hard disks by Alder and Wainwright [1], who found the presence o f the long-time tail in the velocity autoeorrelation function beyond the expected exponential decay on the scale o f the mean collision time, despite many efforts [ 2 - 6 ] , no qualitatively new results have been obtained to date by molecular dynamics. Only recently, by lattice-gas cellular automata [ 7 - 9 ] , Frenkel and coworkers have been able to identify the detailed behavior o f long-time tails. * Corresponding author. Tel.: + 39.59.586059; Fax: +39.59.367 488; e-mail: [email protected]. 0378-4371/97/$17.00 Copyright (~) 1997 Elsevier Science B.V. All rights reserved PII S 0 3 7 8 - 4 3 7 1 ( 9 7 ) 0 0 1 50-7
M. Ferrario et al./Physica A 240 (1997) 268-276
269
On the theoretical side an approach, based on mode-coupling theory, predicts [10], in two dimensions, a power-law decay of the form (v(O)v(t)) o( [(D + v)t] -1 ,
(1)
where D and v are assumed to be constant, and are, respectively, the self-diffusion coefficient and the kinematic viscosity, related to the shear viscosity q by ~/= mpv. The above result, not consistent in two dimensions as it leads to the In t divergence of the, assumed constant, coefficients D and v, has been improved by carrying out a selfconsistent calculation. The self-consistent approach, using time-dependent coefficients [11,12] D(t) and v(t), yields the asymptotic behavior [13]
(v(O)v(t)) cx [(D(t) + v(t))t] -i ~ (t Ix/~-~n t) - I ,
(2)
which gives D(t) cx lx/~7. A similar behavior holds for the stress autocorrelation function, and leads to a slow divergence in two-dimension also for the viscosity coefficient 11. The aim of the present study is to calculate the velocity and the stress autocorrelation functions of a fluid system with a number of particles sufficiently large for the asymptotic regime to be reached and with as much statistics as possible to make an accurate test of these predictions even at high density.
2. Details of the computation We have chosen to study a 2D s o . d i s k fluid, whose particles interact by a pair potential of the form
u(r) = 4~(a/r) 12
(3)
at two state points characterized by the same temperature, kBT/e = 1 and a mass density of pa2= 0.69 and po'2= 0.96, respectively, the latter close to the solid transition [14]. Extensive MD-simulations were performed with a cut-off radius r~. = 1.5 and a time step 6t = 0.005, both values are in reduced m, a, e units, with a corresponding time-unit to = a(m/e) 1/2. The most accurate calculations were performed for a system of 65 536 (64 K) particles and in order to explore system size dependence also systems of 36864 (36K) and 16384 (16K) were investigated. For smaller systems boundary condition effects, related to the recurrence time, interfere severely with the slow, long-time, power-law decay. According to the value of the velocity of sound c = 9.08 calculated by Gravina et al. [14], for the high-density state point the recurrence time is zL =L/c = 28.7 for the 6 4 K particle system, zL = 21.6 for the 3 6 K particle system and zL = 14.4 for the 16 K particle system. At the lower density we expect these times to be even longer, for the combined effects of the increase of the box length and the decrease of the density. The equations of motion were integrated using the velocity-Verier algorithm with on the fly calculation of the autocorrelation functions of the particle velocity Vi(t) and
270
M. Ferrario et al. I Physica A 240 (1997) 2 6 8 ~ 7 6
of the shear stress ~xy = _ _
....
(4)
i:l
j~i
At the same time the mean-square displacement was also directly computed
Ar2(t) =
1
N
~ Z
<[ri(t) -- ri(O)]2> .
(5)
i=1
A large statistical sample was produced using a massively parallel computer, a Cray T3D with 128 processor, with each node computing an independent system trajectory. Starting from a unique initial condition each system was carefully equilibrated by first running it at higher temperatures, with a different, randomly chosen, value of the temperature for each node and subsequently, each system was thermalized applying a sequence of velocity randomizations, using an independent random number generator on each node. During equilibration, trajectories were evoluted according to Nos6-Hoover constanttemperature equations of motion [15]: /'i = P__L,
P~
Pi = Fi -- --~'P',
N
P~ // = -Q-'
t5, = Z
mi
P~ -- -
2NksT,
(6)
i= 1 mi
using a modified velocity-Verlet reversible scheme of integration derived using the Trotter-Suzuki factorization of the Liouville operator [16]. Using a five-term decomposition of the Liouville operator iL, with Ui(z) = exp{iLiz), the approximate discrete time propagator G(t) G(T)-.~ U5(2)U4(2)U3(2)U2(2)Ul(T)U2(2)U3(2)U4(2)U
5
(2)
(7)
is constructed with the following break-up
iL = iLl + iL2 + iL3 + iL4 + iL5 =
F,.N+
.
.
.
.
+ p~. @.
(8)
The resulting algorithm has the following structure, for z = At:
ri(At)
= ri(O) + At ~i(At/2),
r~(At) = {~(At/2) + ~--tmiF~[{r(At)}])exp { - 2 ~(At/2)), q(At) = q(O) + Ate(At/2), At
O(At) = ?7(At~2) ÷ 5-QF~[(~(3t))],
(9)
M. Ferrarioet al./PhysicaA 240 (1997) 268-276 where one first evaluates the half-step values
~(At/2) :
271
¢I(At/2) and ~i(At/2):
r~(0) + ~A--~-t F,[{/'(0)}],
~i(At/2):fi(O)exp {-A---~~(At/2)} + 2A--~tmiFi[{r(O)}]
(10)
and the forces in the above equations are given by Fi[{r(t)}] -
c~U[{r}] , ~ri {r(t)}
F,[{f(t)}] = p~[{/'(t)}] .
(11)
Standard microcanical equations of motion were integrated thereafter and the choice of the velocity-Verlet scheme (or of the modified algorithm used for Nos~-Hoover dynamics) insured the minimum memory usage, since only one 2D position and one 2D velocity per particle needs to be kept in memory at any time. This allowed us to fit large systems on the limited RAM of each node. In fact, the size of the system was primarily limited by memory considerations related to the calculation of the autocorrelation functions along the run, rather than by the computation of the dynamical trajectories. Although systems of more than 100 thousand particles were tested, most calculations were for a system with 16K, 36K and 64K, where K stands for 1024 multiples, this choice being dictated by architectural I/O considerations. Reported results for the 65 536 particle system correspond to a global statistical average over almost 18 millions steps for the low density and more than 30 millions for the high-density point.
3. Results The autocorrelation functions have been calculated for the 64K particles system up to a time t = 40.96 for both state points. At the lower density the velocity autocorrelation function is positive definite in all the time region considered here. It exhibits a shorttime exponential decay, with also a relatively small minimum at t = 0.5 followed by the long-time t -1 regime which approximately starts at time t = 1 and remains very evident all along up to t ~ 4 1 as it is shown in Fig. 1. At the higher density, see Fig. 2, the velocity autocorrelation function has a much more oscillatory behavior at short times and a number of negative-valued minima are present. At around time t = 1 the function becomes positive again and start decaying slowly with a clear t - I power law. For this density the long-time tail is one order of magnitude smaller than for the previous case. Also around time t = 30, that is about the estimated recurrence time, there is a clear evidence of possible interference effects due to periodic boundary conditions. A more careful analysis of the results for the lower density shows some evidence for a somewhat faster than t - t behavior for the larger times. In fact, a least-squares fit of the long-time tail, shown in Fig. 3 gives an exponent o f - 1 within 10 -3, more
272
M. Ferrario et al./ Physica A 240 (1997) 268-276
ill 0.5
0
0
1
2
time
3
4
Fig. 1. Velocity autocorrelation function, p = 0.69. In the inset the long-time behavior is shown using a logarithm scale for both axis. 2
0
0.2
0.4 time 0.6
0.8
Fig. 2. Velocity autocorrelation function, p = 0.96. In the inset the long-time behavior is shown using a logarithm scale for both axis. precisely the exponent is 2 = - 1 . 0 0 0 2 , when the time region between t = 2 and t = 12 is taken into account, while a significantly larger deviation, 2 = - 1.05, is measured when using the far region, from t = 20 up to t = 40.96 (for the entire region a 2% deviation is obtained). Even with the high accuracy and extended range of the present data we feel that it is not possible to discriminate between a genuine faster than t-J behavior and finite-size effects. However, when in Fig. 4 we make a comparison of the velocity autocorrelation functions obtained with different system size for the low-density system, the results practically coincide in the overlapping time region. Notice that the results for the two
273
M. Ferrario et al./Physica A 240 (1997) 268-276
0.004 0.03 0.003
0.01 0.008
0.002
0.006 I
I
2
3
I
I
I
I
I
I ]
20
4 5 678910
30
40
Fig. 3. Power-law fit to the velocity autocorrelation function. The two used time region are explicitly shown, 1.2 0.1
0.8 O.Ol
~
6
K
0.6 0.4
0.001 1
10
0.2
0
l
2 time 3
4
5
Fig. 4. System size dependence for the normalized velocity autocorrelation function, p = 0.69. Curves have been appropriately shifted for a clearer presentation. smaller systems, 16K and 3 6 K particles, are averages performed on a single, long, microcanonical trajectory produced in several weeks on a dedicated workstation. The direct calculation of the mean-square displacements gives results which agree with the slow decay of the velocity autocorrelation function, that is the mean-square displacement grows faster than linearly with time, indeed proportional to In t, as it is clearly shown in Fig. 5 for the system at the lower density. The small deviation, expected according to the previous discussion, is not noticeable on this scale. The situation for the shear stress autocorrelation function is less conclusive. The collective nature of the observable and the overwhelming importance of the potential contribution makes it difficult to obtain high-precision data for the long-time tail, which is expected to be of kinetic origin, and it is not possible to discern a region
274
M. Ferrario et al. IPhysica A 240 (1997) 268-276
0.8 p=.69 0.6
0.4
0.2
0.01
0.1
time 1
10
Fig. 5. Ratio of the mean-square displacement over the time Ar2(t)/t.
8# 105
6 105
104
4 105
2111
103
,
0
,
,
,
I
1
,,,I
,
,
,~-~ 1
,
0.01 ,
0.1 ,
I
I
2
3
time
,
,
,
,
I
,
,
,
,
4
Fig. 6. Shear-stress autocorrelation function, p = 0.69. In the inset a logarithm scale is used for both axis, where a clear power-law decay emerges from the statistical noise. It has to be noticed, however, that, while in the velocity case the magnitude o f the long-time tail in the normalized autocorrelation function is some 3% in the low-density case and 0.2% at the higher density, for the stress autocorrelation functions we do not see any such feature at long times within a statistical noise o f 0.4% at lower density and 0.2% at higher density. The autocorrelation functions for the two systems are shown in Figs. 6 and 7. W h e n the shear-stress autocorrelation functions are integrated over time a welldefined plateau is obtained from which we can extract the G r e e n - K u b o values for the shear viscosity r / = 1.05 + 0.10 at p = 0 . 6 9 and r / = 5.91 4-0.04 at p = 0 . 9 6 , the latter in agreement with the result o f Gravina et al. [14].
M. Ferrario et al./Physica A 240 (1997) 268-276
275
2 106 106
lOS
~k~.
i io 6
0 0
1
2
time
3
4
Fig. 7. Shear-stress autocorrelation function, p = 0.96. In the inset a logarithm scale is used for both axis. 4. Conclusions
Our molecular dynamics calculations are consistent with a high degree of precision with the expected t -1 behavior for the long-time tail of the velocity-autocorrelation function of a two-dimensional fluid. Although some evidence is given for a somewhat faster decay the problem of the nature and the origin of this correction cannot be assessed with the present data. On the contrary, despite the effort, it has not been possible to establish the presence of the long-time tail for the stress-autocorrelation function. There are various possible reasons for that: (i) lack of statistics; (ii) the dominance of the potential contribution to the shear stress, orders of magnitude larger than the kinetic one; (iii) destructive interference of the vortices associated with the individual particle velocity on the overall time memory of the collective observable.
Acknowledgements
We wish to acknowledge enlightening discussions with Daan Frenkel and a critical reading of the manuscript by Ray Kapral. Computer time on the Cray-T3D at the National Supercomputing Center, CINECA, was made available through an initial grant allocation by Cineca and subsequently within the INFM parallel computing initiative.
References
[1] B.J. Alder, T.E. Wainwright, Phys. Rev. A 1 (1970) 18. [2] W.W. Wood, in: E.G.D. Cohen (Ed.), Proceedings of the Summer School on Fundamental Problems in Statistical Mechanics III, North-Holland, Amsterdam, 1975. [3] W.E. Alley, B.J. Alder, Phys. Rev. Lett. 43 (1979) 653. [4] D. Levesque, W.T. Ashurst, Phys. Rev. Lett. 33 (1974) 277.
276 [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16]
M. Ferrario et al./Physica A 240 (1997) 268-276
J.J. Erpenbeck, W.W. Wood, Phys. Rev. A 26 (1982) 1648. J.J. Erpenbeck, W.W. Wood, Phys. Rev. A 32 (1985) 412. M.A. van der Hoef, D. Frenkel, Phys. Rev. A 41 (1990) 4277. M.A. van der Hoef, D. Frenkel, Phys. Rev. Lett. 66 (1991) 1591. M.A. van der Hoef, D. Frenkel, A.J. Ladd, Phys. Rev. Lett. 67 (1991) 3459. For a review, see: Y. Pomeau, P. Rrsibois, Phys. Rep. 19 (1975) 63; see also: M.H. Ernst, B. Cichocki, J.R. Dorfman, J. Sharma, H. van Beijeren, J. Stat. Phys. 18 (1978) 237. B.J. Alder, T.E. Wainwright, D. Gass, Phys. Rev. A 4 (1971) 233. R. Zwanzig, Physica A 118 (1983) 427. M.A. van der Hoef, Simulation Studies of Diffusion in Lattice-Gas Fluids and Colloids, Ph.D. thesis, Utrecht, 1992. D. Gravina, G. Ciccotti, B.L. Holian, Phys. Rev. E 52 (1995) 6123. W.G. Hoover, Phys. Rev. A 31 (1985) 1695. M. Tuckermann, B.J. Berne, G.J. Martyna, J. Chem. Phys. 97 (1992) 1990.