Oscillations in dynamic correlations of binary liquids

ELSEVIER

Physica A 243 (1997) 229-242

Oscillations in dynamic correlations of binary liquids J.C. Lee* Department of Physics and Program in Scientific Computing, University of Southern Mississippi, Hattiesburg, MS 39406-5046, USA Received 2 April 1997

Abstract

Molecular dynamics simulations are performed for a binary liquid mixture (~ +/3) to measure, among others, the collective velocity auto-correlation function Z,~(t) and the time correlation function of the concentration fluctuations F~(t). Both ~-~ and /3-/3 pairs of molecules interact through the same Lenard-Jones potential, but ~-/3 pairs interact only through the repulsive part of the potential. From the results for the relaxation time of Fro(t), we estimate the mutual diffusion coefficients. All measurements probing the individual molecular motions indicate that the molecules execute randomly diffusing motions both at T > Tc and T < Tc. The measurements probing their collective motions indicate otherwise. At T > Tc, the results are as expected from those of the individual motions. At T < Tc, the results indicate cooperative behavior which is not evident in the individual motions. Most importantly, we observe an oscillatory mode of collective motions which is analogous to the optical phonons in two-component crystals.

PACS: 05.40. +j; 61.20.Lc; 66.10. - x

Molecular dynamics simulations contributed significantly to our understanding o f the dynamic correlations in liquids and liquid mixtures. For binary liquid mixtures, previous simulations [ 1 - 6 ] were devoted to compute the mutual diffusion coefficient Dm, using the velocity auto-correlation function formalism o f Green and Kubo or the equivalent Einstein formalism [ 7 - 9 ] , and to compare it with the self-diffusion coefficient D. The velocity-time correlation function provides a useful piece of information about molecular motions, but it cannot be measured in the laboratory. The time correlation function of the concentration fluctuations can be measured by performing a dynamic light scattering experiment [8], but to the best of our knowledge there has been no effort to compute it to date. In this paper, we will compute both the collective-velocity * E-mail: [email protected]. 0378-4371/97/$17.00 Copyright (~ 1997 Elsevier Science B.V. All rights reserved

PHS0378-4371(97)00281-1

230

J.C. Lee/Physica A 243 (1997) 229 242

auto-correlation function Zm(t) and the concentration-fluctuation auto-correlation function Fro(t). The system we simulate is a binary liquid mixture which exists in a homogenous mixed phase at high temperatures and in two separated phases at low temperatures. At temperatures in the one-phase region, our results for Fro(t) are as predicted by the hydrodynamic theory and therefore it is possible to estimate the mutual diffusion coefficient. In order to determine D m accurately one has to probe the hydrodynamic limit of long wavelength which requires a large system. Although our system is larger than those explored previously, it is still relatively small, and therefore it is not possible for us to speak of the true hydrodynamic limit. However, the results for long wavelengths show that the wave-vector-dependent diffusion coefficient Din(q) becomes nearly independent of q, which we consider sufficient for the purpose of finding out the general trend of Dm and D. In the two-phase region, on the other hand, both functions exhibit relaxation patterns which suggest cooperative behavior. We observe a quite convincing evidence of a normal mode which is analogous to the optical phonons in crystals with two different atoms per primitive basis. The system to be simulated consists of N = 4096 Lenard-Jones molecules in a box of L × L × L subject to the usual periodic boundary condition. Half of the molecules are species ~ and the other half are species ft. Both ~-~ and fi-fl pairs of molecules interact with each other through ~b(r)= 4~{(tr/r) 12 -(tr/r)6}, while all ct-fl pairs of molecules interact through q~(r)= 4e(a/r) 12. Both species are assumed to have the same molecular weight. To contrast the two species in the concentration fluctuations, we assign ÷ 1 to the dielectric constant of all molecules of species a, and - 1 to that of all molecules of species ft. The Verlet leap-frog algorithm is adopted with 6t = 0.004. The time variable t is given in units of (mtr2/48e) 1/2, and all length in units of a. Choosing L = 17.03, the density is n =0.83. What we will call temperature T is actually kT/e= (16/N)~i(u2i). The total momentum remains zero. For each T, measurements are preceded by a long thermalization process to ensure that the system is in equilibrium. We will compute three pairs of dynamic functions which reflect how molecules are moving around in the mixture. In each pair, the first function attempts to describe the effect of the individual motions on the pertinent dynamics while the second one attempts to describe that of the collective (mutual) motions [7-9]. The first functions may be computed quite accurately because one can take an average over the results for N different tagged molecules. For the second functions it is not possible to take such a large number of samples, and therefore it is much more difficult to determine them accurately. It is even more so at T < Tc because the time scales necessary to describe the second functions are larger than those at T > Tc by almost an order of magnitude, which demands a great deal of CPU and computer memory. The first pair are

(1)

J.C. Lee I Physica A 243 (1997) 229-242

231

and (2) The first function (r2)~ is the mean squared displacement of the individual molecules of species ~ while the second function (R2)~ is the mean-squared displacement of what we will call conveniently the center of mass of molecules of species a. We did not divide Eq. (2) by N2/4 as one should for the true center of mass. In this way, (r2)~ is the self-part of (R2)~, and the difference between the two functions (R2(t))~ - (r2(t))~ = Z

({r~i(t) - r~i(0)) • {r~j(t) -

r~j(O)})/(N/2)

(3)

i¢-j gives a measure of the cross-correlations among the displacements of different molecules. Eqs. ( 1 ) - ( 3 ) are computed separately for molecules of species ft. As expected, the results for the two species turn out to be the same. Thus, the two results are averaged for the final results which we will present without subscripts. Because (r 2) = 6Dt in the asymptotic regime, the self-diffusion coefficient D can be measured easily from the results for (r2). In a similar fashion, one can, in principle, measure the mutual diffusion coeffÉcient Dm from the measurement for (R2), but in order to determine Dm in this way one also has to compute the thermodynamic factor Q = (x~/kT)(Sp~/Ox~)r,p, where x and # are the mole fraction and the chemical potential, respectively [11]. Using the Kirkwood-Buff theory [12], one can reduce Q to a set of very simple integrals over various radial distribution functions [2,3]. Unfortunately, these integrals represent the fluctuations of various products of particle numbers and therefore call for a simulation with an open system which the present system is not. For this reason, we make no effort to compute Q. The second pair consist of two velocity auto-correlation functions defined by

Zs(t ) = ~'~i (U~i(t) " uai(O)) ~i (Hcti(O)-'~)

(4)

and

( U ( t ) . U(O)) Zm(t)= (U(O)" U(O)) '

(5)

U(t) = Z u~i(t) - Z u#i(t).

(6)

where

i

i

Here again Z~ is computed separately for the molecules of species fl and then the two results are averaged. Due to the momentum conservation law, U ( t ) = 2 ~ i u~i(t), and therefore Z~ is the self-part of Zm.

J.C Lee/Physica A 243 (1997) 229-242

232

The third pair describe how the molecular motions affect the concentration fluctuations (intermediate scattering functions):

(Ps(q, t)ps(-q, 0))

Fs(q, t) = (ps(q,O)p~(-q,O))

(7)

(p( q, t ) p ( - q , 0)) Fro(q, t) = (p(q, 0 ) p ( - q , 0)) '

(8)

and

where Ps( q, t) = exp(iq, ri ) and p( q, t) = ~ i exp(iq, r~i) - ~ i exp(iq, rl3i). The negative sign in the second term of p is due to the dielectric constant - 1 that we assign to the molecules of species ft. The first function G(q, t) probes the concentration fluctuations due to the individual molecular motions. It is given by [7-10] G(q, t) = e x p ( - F ( q ) t ) ,

(9)

where the relaxation rate F(q) is related to the self-diffusion coefficient D by

F(q)=q2D.

(10)

The second function Fm(q, t), on the other hand, probes the concentration fluctuations due to the collective motions of molecules. Assuming that the thermal diffusion coefficient is much larger than the mutual diffusion coefficient, the hydrodynamic theory [ 13,7,10] predicts that

Fro(q, t) = e x p ( - F ( q ) t ) ,

(11 )

where the relaxation rate F(q) is now related to the mutual diffusion coefficient Dm by

F(q) = q2D m .

(12)

Eqs. (11) and (12) assume that the wave vector q is small enough to probe the hydrodynamic limit. When q is not small enough to do so, one can only speak of a q-dependent effective mutual diffusion coefficient Din(q). Writing the wave vectors in the form of q = (2niL)k, we choose k = 2, 3, 4, 5, 9. (Because some of these choices for k were made at different stages of the computation, not all quantities were computed for all of these wave vectors.) The corresponding q will be referred to as q2,q3, etc. The measurements are taken with each wave vector directed along the x - , y - , and z - a x i s , and then the three results are averaged. For each q, the probed length scale is given by 2 = L/k. Two comments are in order concerning the periodic boundary condition. When a molecule leaves the box in which it was initially placed (the primary box), the periodic boundary condition brings the image of the molecule back into the box on the opposite side. This action should not make a false contribution to the dynamics that we measure. For that reason, when we compute (r 2) and (R2), we follow the molecules as they leave, or as they retum to, the primary box. The only effect of the periodic boundary condition is that when the molecules enter the next box, they find themselves

J.C. Lee/Physica A 243 (1997) 229-242

O

• O0

•

•

OO

•

• 0 0

•

o

•

0

•

0

o 0

O0

0 0

O

0 0

•

•

0

o

0 •

•

O0

• •

0

•

• •

•

0

•

O0

•

•

Oo

0

•

• o •

0 O0

0

0 0

0

• •

0

•

• •

• 0 0 0

•

o 0

0

0

•

•

o

0 0

e •

0

0

•

•

D

0 0 0 0 o 0 o

•

•

©

•

e0

0

•

0

• •

0

0

o

0 0

•

0 •° 0

0

0

•

O0

0 i)

0 •

•

0 •

•

•

o

•

•

0

o0

0

•

o

0

0

•

0

o

0 0

0

0

0

•

•

•

•

0

•

0 •

0

Q

o

•

•

0 0

233

0 0

0

0

0 o

•

•

•

•

0

0

• ° ° °e

°

•

•

® •

•

0

• 0

•

0

0

D

O0

0 9

Fig. 1. Snap shots of molecules with x-coordinate anywhere between 8.6 and 9.4. T = 28.12.

in the same environment as their image molecules do in the primary box. The periodic boundary condition poses a much more serious pitfall in the computation for the relaxation of the concentration fluctuations. When there is an excessive local concentration, the excess can only relax by transferring it over the distance of the probing wavelength, which is why Eqs. (11) and (12) predict that the relaxation time should increase with the increasing wavelength. Because the periodic boundary condition instantaneously transfers molecules over the distance of L, it is necessary that the probing wavelength is shorter than L. That is why k = 1 was avoided in the above. First we present the results for T=28.12. Shown in Fig. 1 are the molecules that have x-coordinate anywhere between 8.6 and 9.4. It is clear that T > Tc and the system is in the mixed homogeneous phase. Fig. 2 shows the results for (r 2) and (R2). To obtain these results, our average was taken ever approximately 80 different initial times where two successive initial times are separated by t -- 12. This was repeated four times, each time starting from a different initial configuration for a different thermalization. The figure shows that (R 2) is approximately half of (r2), and therefore the displacements of different molecules tend to be negatively cross correlated. The negative sign indicates the absence of cooperative behavior in the individual motions. The slope of the (r 2) curve yields D = 0.204. Fig. 3 shows the results for the velocity auto-correlation functions. The two autocorrelation functions relax following the same pattern, which means that the motions of the individual molecules and those of the two centers of mass are essentially the same except in magnitude. This seemingly reasonable expectation is not met at T < Tc, as will be shown later.

J.C. LeeIPhysica A 243 (1997) 229-242

234

0

R2

000[30 rz

0

d-

[] O nOD O

o_ ,I

[] [] D

[]

[]

[] []

[]

[]

[]

00000

C3

ODD oo ODDD

O_

u~

,...0000~ 0 Q'-' 0oO OO

DE~OO oo©O000

0

°0.0

5.0

10.0

Fig. 2. The m e a n - s q u a r e displacements at T = 2 8 . 1 2 .

o.- 0

cx:o~ z~(O rToooD Z.(O

~o

0

6-D

0 ~Jd N

0 [] 0 0

O_

_

_

d (N

IO.O

'

5.10

'

'

[] D

' 1 0 ~.0

t

Fig. 3. The velocity auto-correlation functions at T =28.12.

J.C. Lee/Physica A 243 (1997) 229-242

235

\ \ \ \\ \

EF LL - ......

\

_

J

0 .lO

~

I

I

1

5.0

*

q=q2 c:qa

c=q9

_

I

l

I

i

10.0

f Fig. 4. The auto-correlation function Fs(q, t) at T = 28.12. The values of the three wave vectors are given in the text.

Fig. 4 shows the results for F~.(q,t). The results for all wave vectors follow Eqs. (9) and (10). The relaxation rates are F(q2)=O.11, F(q3)=0.25, F(q5)=0.68, and F ( q 9 ) = 1.88. Eq. (10) gives D = 0 . 2 0 , in agreement with the result from (r2). Fig. 5 shows the results for Fro(q, t). The exponential decay pattern is as predicted by Eq. (11 ). We find that F(q2) = 0.07, F(q3 ) = O.16, F(qs) = 0.51 and F(q9) = 1.85. Eq. (12) then gives: D,~(q2) = 0.13, Din(q3) = O.13,Dm(qs) = 0.15, and Dm(qg) = 0.17. The numerical uncertainty does not allow more than two significant digits. Within this limitation, we estimate that Dm= 0.13, or more safely that Dm is less than 0.13. This is significantly less than D, and it is quite clear that the relaxation does not follow the Hartley-Crank approximation [14] which states that Dm----xIDI + x 2 D 2 where x and D represent the molar fraction and the self-diffusion coefficient, respectively. This is due to the repulsive inter-species interactions. When random fluctuations create excessive concentration of either species at any place, the repulsive inter-species interactions make the excess to relax rather slowly even in this one-phase region of the phase diagram. Next we compute the same set of functions for T = 6.25. Fig. 6 shows that the system is now in the two-phase region of its phase diagram. Fig. 7 shows the results for (r 2) and (R2). It is notable that (R 2) is now larger than (r 2) at all times. This means that the displacements of different molecules tend to be positively cross correlated and therefore the center of mass moves faster than the individual molecules do. Fig. 8 shows the results for the velocity auto-correlation functions. It is notable that Z,, relaxes much more slowly than Z~, suggesting collective behavior. The slow

J.C. Lee/Physica A 243 (1997) 229 242

236

\

\

""'''-.

\ \

" ' ' ' - .

\

6-

\ \ \

I,

\

E

\ \ \

0

\ \

\ \ \

I

- " "

--- -

0.'0

'

'

5 .lO

,

q=q2 q :qs q=% q=q9

, 10~,0

'

f Fro(q,t)

Fig. 5. T h e auto-correlation function

•

•

©

•

•

0

0

0

0

O

0

•

0 O

0

•

0

0

0

•

O

•

Q

•

•

•

® •

•

O

•

•

•

®

• Q

•

@

•

•

•

g,

Q

@

•

•

• •

• •

• •

•

•

• •

•

• -

-

•

• 0

•

•

•

9

•

•

• •

•

@

• •

O

•

0

• •

o

0

0 •

0

©

0

•

o •

0 'Z

• •

0

0

0

0

0 0 0 0 0

O0 0

0

0

0

•

0

0 O0 0

o

• •

•

~

•

O

0 0 0

©

•

•

0 0

~

•

•

0 0

0

Q o°

o •

0 0

0

0

• qD

0 0

•

•

•

0

0

•0

• 0

0

OOo

0

0

•

0

0

0

•

•

at T = 28.12.

•

•

•

gB

Fig. 6. Snap shots o f m o l e c u l e s w i t h x - c o o r d i n a t e a n y w h e r e b e t w e e n 8.6 a n d 9.4. T = 6.25.

J,C. Lee/Physica A 243 (1997) 229-242

237

0

60

O o

O00CO

R2

o.

00000

r2

J 0 ~

o

0

oo

O_

0

0 0 ©

C3

o

/

d

-

oo°

/

oj

cm

do.T-

'

~c~

20',0

'

4o'.o

'

60'.0

'

Fig, 7. The mean-square displacements at T = 6.25. o-

0

z,(.0

.....

Z,.O)

ozooo 0

+-~.

-~,0

e~

oO_

6

0-4

6~0.0

'

50~,0

'

, 100.0

,

-m 150.0

f Fig, 8. The velocity auto-correlation functions at T = 6 2 5 .

J.C. Lee/Physica A 243 (1997) 229-242

238

relaxation is superimposed by an oscillating pattern. The oscillation is remarkably periodic and is present regardless of the initial configuration before the therrnalization. One may suspect that this is an artifact due to the periodic boundary condition. The period of oscillation is approximately 11 time units. During this period the center of mass moves approximately 2 length units (2a). This is too short to have any thing to do with the periodic boundary condition. To explain what the oscillation may then mean, consider the external field F(t) that interacts with the molecules as shown by the Hamiltonian

~'=

(13)

- A(t). F(t),

where A ( t ) - - ~ i r ~ i written as

~'~irlji. When this force oscillates, its effect on U(t) may be

(A U(t)) = ZuA(~o)Foexp(-io~t), where the dynamic response function

(14)

lUA(~O) is given by

ZuA(og) = f l / ( U ( t ) . U(0)) exp(i~ot) dt.

( 15)

o ,

0

Because F(t) is coupled to A(t), the dispersion and dissipation is determined by the dynamic susceptibility of the system to the force field, ZAA, not ZtJA(Og). The two functions are, however, related to each other by ZUA(Og)=

flCuA(t ----0) -- i~OX,~A(CO) ,

(16)

where CAA(t)= (A(t). A(O)). Thus, the real part of ZtJA(Og) gives the dissipation and the imaginary part the dispersion. The general pattern of ZuA(~o) is shown in Fig. 9. The peak located approximately at 09o = 0.6 accounts for the oscillating pattern of Zm(t), and suggests the existence of a normal mode of collective motion. It is tempting to regard this as the neutral-particle analogue of the plasma oscillation of the charged systems, but the present system lacks two essential elements of the charged systems necessary for the charge oscillation. The interaction is not long ranged and the system has no tendency to neutralize the local "charge" distribution, i.e., there is no tendency for the two components to mix. It is pertinent to note that the oscillation of Zm(t) disappears as the temperature is raised toward Tc and therefore the normal mode exists only at temperatures in the two-phase region. This leads us to suspect that the oscillation arises from the interface. As the interface tries to smooth out itself, the two centers of mass move, but the molecular inertia overcorrects the shape and the process repeats. It is a liquid analogue of the optical phonons. Fig. 10 shows the results for F~. The relaxation is again consistent with Eqs. (9) and (10), and we find D=0.055. No noticeable change is observed which may be attributed to the phase separation. By contrast, as Fig. 11 shows, the results for Fm are quite different from those we observed at T > Tc. The results with q = q2 deserve a particular attention. Here the concentration fluctuation is probed in the length scale

£C. LeelPhysica A 243 (1997) 229-242

239

00

Xi )

c5tO (D

O

04 O

(a)

O f

o 0.0

--

t

1

t

1

~

I

u

I

0.5

1.0

0.5 CO

1.0

OO

c;-

/¢

X(~0) tO t:D

,5c5O_

o 0.0 (b)

Fig. 9. (a) The real part of ZUA(OO).(b) The imaginary part of ZuA(co). comparable to the domain size, and therefore the results can provide an information about the collective domain motions and they do. There is an oscillating pattern in the q2 curve. Although there is a considerable degree of numerical uncertainty, we observe this oscillation regardless of the initial condition before the thermalization. The period of oscillation is approximately 11 time units, which makes us to believe that this is the same mode of collective motion observed earlier in the Zm data. The oscillation occurs somewhat more prominently at T = 5.08, as shown in Fig. 12. To confirm that all of the above cooperative types of molecular behavior at T < Tc are a result of the phase separation, we momentarily replace the soft-sphere repulsive interactions (between molecules of different species) with the full Lenard-Jones potential, and repeat the above computation. The temperature is chosen to be the same as above, i.e., T = 6.25. In this case there is no phase separation. The replacement has no noticeable effect on (r2), Fs, and Zs. However, all the dynamic functions pertinent to collective motions do change. The (R2) curve now appears below the (r 2) curve, and Zm becomes identical to Zs, and there is no oscillation in Zm, nor in Fro(q, t). The results for Fr~(q,t) shows that Dm is now larger than D, which confirms our earlier interpretation for the result, Dm < D at T > Tc. To summarize, at temperatures in the one-phase region we have shown that it is possible to probe the concentration fluctuation near the hydrodynamic limit and estimate

J.C. LeelPhysica A 243 (1997) 229~42

240

\\'.. \

"-.

1

\

\

"'-

\

"'-.

\

k

"-.

\

\

\

"-.

\ %

\

". -.

\

"-.

%

\

"\

\

-, "''-

1

\

.

\

\

\

1

%,

\

\ %,

1

O9

\

L

~I

- ...... ---0.'0

'

--

\

\

\

q=q2 q=q~ q=q4 q=qs q=q9

I 0'. 0

'

x,

,,, \

,,

20.0

3C,0

'

40r.O

f Fig. 10. The auto-correlation function ~ ( q , t ) at T = 6 . 2 5 .

"~

"'''''-4.

,%

.+.-

U" L

/

"" --.

E 0

\

- ......

\

-

1

-

----

-

q=q2 q=qs q=q4 q=q9 q=q5

Sc5:

0/0

10.0,

,

20.0I

,

3010

'

t Fig. 11. The auto-correlation function Fro(q,t) at T = 6.25.

4o'.o

241

J.C. Lee/Physica A 243 (1997) 229-242

..+.-

Oi,

E <5-

-...... - - -

q=q2 q=q~ q=q4

0.

0.'0

10'.0

20'.0

30'.0

40'.0

Fig. 12. The auto-correlation function Fm(q,t) at T = 5.08. the mutual diffusion coefficient. At temperatures in the two-phase region, we have observed a reasonably convincing evidence of an oscillating mode of collective motion. The intermolecular forces in most of the phase-separating binary liquid mixtures are more complex than assumed in the present model system. Nevertheless, the oscillatory mode should be present at temperatures in their two-phase region. Our original purpose was to investigate the dynamic correlations in the presence of a quenched random disorder. However, it was necessary for us first to understand the dynamic patterns of the pure system before we could speak of the effect of the disorder. It is hoped that we will be able to report in the future how the dynamic correlations reported here are altered in the presence of the disorder. The author wishes to thank G.S. Grest and R. L. Gibbs for their help while he was writing the MD code. Helpful discussions with D. McCain, V. Yamakov, G. Rayborn, and D. Stauffer are gratefully acknowledged. The computation was performed on a CRAY/YMP at the Mississippi Center for Super Computing Research. This research was supported by the Donors of the Petroleum Research Funds administered by the American Chemical Society.

References [1] G. Jacucci, I.R. McDonald, Physica A 80 (1975) 607. [2] D.L. Jolly, R.J. Bearmann, Mol. Phys. 41 (1980) 137. [3] M. Schon, C. Hoheisel, Mol. Phys. 52 (1984) 33.

242 [4] [5] [6] [7] [8] [9] [10] [11]

J.C Lee/Physica A 243 (1997) 229 242

D.M. Heyes, J. Chem. Phys. 96 (3) (1992) 2217. J. Trullas, J.A. Padro, Phys. Rev. E 50 (1994) 1162. Y. Zhou, G.H. Miller, Phys. Rev. E 53 (1995) 1587. see, for example, J. P. Boon, S. Yip, Molecular Hydrodynamics, Dover, New York, 1980. J.P. Hansen, I.R. McDonald, Theory of Simple Liquids, Academic Press, San Diego, 1986. D.A. McQuarrie, Statistical Mechanics, Harper Collins, New York, 1976. B.J. Berne, R. Pecora, Dynamic Light Scattering, Krieger, New York, 1990. See, for example, D.D. Fitts, Nonequilibrium Thermodynamics, McGraw Hill, New York, 1962; H.J.M. Hanley in: H.J.M. Hanley (Ed.), Transport Phenomena in Fluids, Dekker, New York, 1969. [12] J.G. Kirkwood, F.P. Buff, J. Chem. Phys. 19 (1951) 774; see, for a review, H.L. Friedman, A Course in Statistical Mechanics, Prentice-Hall, Englewood Cliffs, NJ, 1985. [13] R. Mountain, J. Deutch, J. Chem. Phys. 50, 1102 (1969); H.N.W. Lekkerkerker, W.G. Laidlaw, Phys. Rev. A 7 (1973) 1332. [14] R.J. Bearman, J. Phys. Chem. 65 (1961) 1961.