THE UN I V E R S I T Y OF M I C H I G A N COLLEGE OF ENGINEERING Department of Nuclear Engineering Technical Report FLUCTUATION ANALYSIS IN SIMPLE FLUIDS A. Z. Akcasu E. Daniels ORA Project 08034 supported by: NATIONAL SCIENCE FOUNDATION GRANT NO. DK-1009 WASHINGTON, D.C. administered through: OFFICE OF RESEARCH ADMINISTRATION ANN ARBOR July L969

This report was prepared as a paper for submission to The Physical Review. ii

TABLE OF CONTENTS Page ABSTRACT 1 1. INTRODUCTION 2 2o THE GENERALIZED LANGEVIN EQUATION 5 3. TRANSVERSE CURRENT CORRELATIONS 8 4o A GENERALIZED HYDRODYNAMIC DESCRIPTION 16 Transverse Equations 29 Longitudinal Equations 30 Conclusions 36 ACKNOWLEDGMENTS 37 REFERENCES 38 FIGURE CAPTIONS 40 iii

ABSTRACT Generalized Langevin Theory is applied. to the analysis of fluctuations in simple fluids. Formulas for the current-current and density-density correlation functions are developed. More importantly, they are developed for a region of frequencies (a - 1013sec 1) and wave vectors (k 108cm-1) which are explored in typical slow neutron scattering measurements. Where applicable, comparisons a d re made ith the results of the numericallculations carried out by Rahman, and good agreement is generally found. 1

1. INTRODUCTION The calculation of correlation functions in classical simple fluids in terms of microscopic quantities plays an important role in statistical mechanics, both for the interpretation of scattering experiments and the evaluation of the frequency and wavelength-dependent transport coefficients. Among these, the density-density (or Van Hove) correlation function has received the most attention because of its direct relation to the differential scattering crosssections. Other correlation functions, such as the transverse current-current correlation functions which are not readily accessible in experiments, have been subject to quantitative investigations only recently following the publication of Rahman's) computer calculations in argon-like liquids. Computer (2) studies of correlations using molecular dynamical calculations provide a stringent test of the validity of the various classical theories introduced in correlation analysis because they only assume a known model interparticle potential, and involve no quantum effects. The classical analysis of correlations is usually based on either a kinetic or hydrodynamic description of fluids. The kinetic description developed extensively by Nelkin and his co-workers has been justified theoretically for (4) dilute gases, and used successfully to interpret-Brillouin scattering from gases (5,6) It has also been extended. to dense fluids and applied. to Rahman's (2) molecular dynamics calculations for liquid argon) with poor quantitative agreeent (8) The hydrodynlamic description of fluids has long been in use in the fluctua2

tion analysis in arbitrary continuous media as a phenomenological theory. In this approach, the conventional hydrodynamic equations are used to describe the linear response of the fluid, and the correlations functions are then related to the linear response by means of the fluctuation-dissipation theorem. (10) (11) Using the formulation described by Landau and Lifschitz, Rytov applied the fluctuation-dissipation theorem to distributed parameter systems and calculated, among others, density-density correlation functions in an arbitrary continuous medium. This hydrodynamic approach has been used to interpret light scattering from liquids successfully by Mountain. A systematic and general hydrodynamic description of fluids for the calculation of correlations and (12) transport coefficients has been developed by Kadonoff and Martin. The (1)4) latter approach has been applied by Chung and Yip to Rahman's calculations of current-current correlations. The objective of this paper is to present a classical analysis of correlations in simple fluids, based on the generalized Langevin equation developed by Zwanzig, and Mori and to interpret quantitatively the current(2) current correlations computed by Rahman for liquid argon. This approach has several appealing features. First by choosing the dynamical variable in the description of the fluid as the microscopic phase density function one obtains ) an exact kinetic equation for the correlation function r(v,v';x,x',t), (19) which reduces by approximation to the kinetic equation derived by Zwanzig (8) and Nelkin. However, by making an alternative choice of the dynamical variables to be microscopic densities in configuration space (e.g., mass, current, and energy densities) one arrives at an exact hydrodynamic description of cor5

relation functions in terms of frequency and wavelength-dependent transport parameters. The choice of the appropriate set of dynamical variables is arbitraryo For any choice of these variables one obtains exact expressions for the correlation functions of the variables in the set. Different levels of approximations can be obtained. for a particular correlation function by adding new variables to the set and. using the same simplifying assumption (e.g. Markoff assumption) in each case. The continued fraction expansion of correlation (17) functions by Mori, for example, can be obtained by using an orthogonal extension of the-set starting from a given dynamical variable. The separation of thermodynamic and transport parameters can be given a geometric interpretation in terms of projections of dynamical variables on appropriate orthogonal axes, and the extension of their definitions to short wavelength where the anisotropies become significant can be made in a systematics.way. In this paper we will use the configuration space (hydrodynamic) description of fluids to investigate the current-current correlation functions, and obtain approximate formulas for transverse and longitudinal current-current correlations in the frequency and wavelength regions encountered in neutron scattering using a -;ov..........,v assumption. The results will then be compared to (2) Rahman's computer data for liquid, argon. We will also obtain the wavelength and frequency dependence of shear and longitudinal viscosities explicitly and discuss the influence of thermal effects as a function of wavelength. The objective of this paper is similar to that by Chung and Yip; however their approach is based on the correlation function formalism developed by Martin and Kadonoff, rather than on the projection operator formalism by Zwanzig and Mori. 4

2. THE GENERALIZED LANGEVIN EQUATION Extending the projection operator technique first introduced by Zwanzig, Mori (6) proves that the equation of motion of a. set of dynamical variables, aj(t), can be written in the form, dt) - i-.a(t) + t cp(t-u).a(u)du = f(t), t > 0 (2.1) The state vector a(t) is defined such that it has no invariant contribution, e.g., a(t) = A(t) - <A(t)>, (2.2) where <...> denotes the thermal average of the vector A(t). Equation (2.1) is the generalized form the Langevin equation (10,20) in the stochastic theory of Brownian motion. The random force vector f(t) is given formally by f(t) = et(-P)iL( ), (a = iLa), (2.3) where L is the classical Liouville operator and P is a, projection operator defined for any arbitrary phase function G(t) by PG(t) = <G(t)a >.<aa t> a. (2.4) t Here, a denotes the row vector which is the hermitian conjugate of a. The t - t <aa > is the inverse of the square matrix [<aiaj>] which is the static correlation matrix. It should be noted that the evolution of the force vector f(t) is determined by the special propagator e, whereas the evolution of the state vector a(t) is given by 5

a(t) = e a. (2.5) It is shown by Mori that <f(t)a> 0, t > 0. (2.6) The square matrices cp(t) (damping function) and Q (frequency matrix) are defined by p( (t) < f(t)(o)>.<aa>-l, (2.7) and = <aa >.<aa>. (2.8) Multiplying (2.1) by a *<aa >-1 from the right, taking the thermal average of the resulting equation and using (2.6) yields dR(t) iQR(t) + t c(t-u)-R(u)du = O, t (2.9) dt o where R(t), the normalized dynamic correlation matrix, is defined by R(t) = <a(t)a'><aa >-1. (2.10) The one-sided Fourier transform of R(t) is obtained from (2.8) as R(iw) = [iw-iO+cp(io)]-, (2.11) where C(ic) = lim f e (i )t (t)dt. (2.12) 6+0 The projection technique enables one to find a closed set of linear equations for the correlation matrix R(t) when the state variables a.(t) are chosen 6

as fluctuations from thermal equilibrium. The theory, although formally exact, serves only to transform the calculation from the direct computation of R(t) to the computation of Q and cp(t). However, the frequency matrix Q is determined from static correlations which are generally much easier to compute than the time-dependent correlation functions, and we may use approximations to compute the damping matrix cp(t)^ In particular we will consider representations in which we can make a Markov approximation on cp(t), viz., = lim cp(ic), (2.13) woO to approximate the transform of the correlation matrix R(iwo) by R(iw) = [ia)-inm+]1. (2.14) We shall follow this formal procedure to calculate the current-current correlation function and its transform by choosing the components of the state vector as the spatial Fourier transforms of the local densities of conserved variables. 7

3. TRANSVERSE CURRENT CORRELATIONS We first calculate the cosine transform of the transverse current correlation function as a simple application of the generalized Langevin equation, and then compare it with Rahman's computations. For this purpose, we choose the components of the state vector as al = Jl(k) (3.la) a2 = 31(k) (3.lb) where J(k) and H1(k) are the mass current density and the stress tensor, respectivelyo They are defined by N J-(k) - mvjexp(ik.x), (3.2a) N r N x.x. ^ *= x -Ax (5.2c) ap dV(R) l-e- ik.R P (k ~ R dR ikR R ap 3(k 2d) dR ik.R In these definitions, x and v denote the position and velocity of the ath particle in the system, and the subscripts i and j refer to the Cartesian components in a coordinate system in which k is parallel to the z-axis. With this choice of variables the static correlation matrix is diagonal, viz., 8

< a <alal> 0 - <a a > = (33) L 0 <a2a2>J since the variables al and a2 are respectively odd and even functions of the particle velocities so that <a2al> = 0. Furthermore, a direct evaluation of the diagonal terms using (3.2) yields a <ala=> p (3-4) <a,2a2> = 44(k) (3- j( [n + n2fd3R g(R) a2V (1-coskZ)]:gR aX2 ka where 5 = (l/kBT)(kB= Boltzmann's constant), V is the volume of the system, and p is the equilibrium density. Here C44(k) is one of the elastic moduli cal(21) culated by Schofield. The frequency matrix S; is calculated using (5.h) and (35,), and noting al = ika2 as follows: 0 1 Q - kL j (3.6) 1544(k)/P0 oJ Since (l-P)al = ik(l-P)a2 = 0, the random force component fl(t) is identically equal to zero and the damping matrix cp has only one non zero element, p(t) =. (3.7a) L0 0p22( where p22(t) = <(1-P)a2 e(1-P)al <a2a2>1 (3-7b) 9

The generalized Langevin equation for the set (al(t),a2(t)) becomes da(t) - ika2(t) = 0 (3.8a) dt da2 - ik a al(t) + t q22(t-u)a2(u)du = f(t). (8b) dt p o 0 The transverse current correlation function is defined by R1(k t) <a (t)a> (9) =kI a (3.9) <alal> Its cosine transform R1(k,w) is the transverse current power spectral density. The latter can be obtained directly from (3.8) by multiplying it by al, taking the thermal averages, using <f2(t)al> 0, and Rl(k,n) = Re[iw+(k2/p )r (k,iw)]'1 (3.10) where T (k,io) may be identified as a k- and. -dependent shear viscosity defined by Po <a2(io)al> T(kiw) =- <a a- (3.11a) <ai(iiw)al> C44(k) (3.11b) iWc+Cp22(k, iW) (2) Rahman() has computed RL(k,o) for various values k and a. We have been able to obtain an exact expression for it in terms of the Laplace transform of the damping function, viz., ~22(k,iw). However, the evaluation of the latter using (3o7b) is as difficult as solving the Liouville equation, although perturbative 10

(18) techniques such as expansions in density or interparticle potential may be used in dilute systems. For sufficiently small frequencies, we can approximate (3510) by replacing p22(k,iw) by its zero frequency limit 22(k,o). This approximation corresponds to a Markov description of the fluid in terms of al(t) and a2(t), in which the convolution integral in (3.8b) is replaced by ot p22(t-u)a2(u) = a2(t) f du p22(u). (3.12a) The frequency range in which the Markov assumption may be expected, to be valid can be estimated by considering the next term in the expansion of p22(k,iw) in powers of (iw): W << cp22(ko)/d 2(ic) iw=o (35l2b) It is clear that the Markov assumption ceases to be valid if 22 (k,o), the leading term in this expansion vanishes for some values of k. In such cases it turns out to be more convenient to go to a more complete description of the fluid by introducing new variables. With these remarks, we obtain the following approximate form for-rR(kc): W (k)k2C44(k)/p R_(k,wo) S - (3513) 2 2 k2C44(k) 2 C0 OD (k; + cW S P where we have introduced ws(k) = lim P22(k,_1) ~ (3.14) 11

Calculation of cs(k) directly from (3.14) is still a formidable task in dense fluids, although it is simpler than calculating c22(k,icU). Therefore, we choose to try to guess its k-dependence by considering the asymptotic behavior of RL(k,c) in the large-k limit, and of s(k,iw) in the small-k limit. First, it is known that lim V (k ) 4() (3.15) k,o w5(o) (5.15(o) k,; +o s is the conventional shear viscosity, s. Hence, the small k limit of cos(k) is given by (o0) = 44(o (3-16) s (Note that C44(0) is G in Zwanzig's notation(22). 00 The large-k behavior of R_(k,w) may be predicted from the transverse current spectrum of an ideal gas, viz., IG _ - m2W/2k2 R1 (k,) - d e- 2/2 (35.17),2k2/ This function has a single maximum at cu = o for all k. On the other hand, (3.13) attains its maximum for a fixed k at a frequency < (k) = kC4 ) (k) (3.18) sm p 2 for o2(k) < 2k2C44(k)/p, and at C = o otherwise. (The function u (k) is often S 0 sm referred to as the "dispersion relation. ")() The dispersion relation (3.18) is expected to approach that of an ideal gas as k a- Xo because for large-k the par12

tides behave as free particles. (This is more apparent in the case of longitudinal current correlation function because it is related directly to the neutron scattering cross section, there the large values of k correspond to large momentum transfer to the scattering medium.) Hence, we require e2 (k) to approach sm zero for large values of k, i.e., 2 2k2 44(k) 0C2(k) (k (3.19) s P 0 Thus, we obtain the asymptotic behavior of ws(k) for small and large k from (3.16) and (3.19). The k-dependence of %s(k) for the intermediate values of k can be obtained by interpolating it between the zero and large k limits by the following formula: 2k C44(k) + L )22(C44(k)-pJ/Pm)/P]. ob (k) + (3.20) o [l+k2/k2] where k is an adjustable parameter whose choice will be discussed presently. It is interesting to compare (3.13) using the above expression for es(k) to the ideal gas result in (3.17) for zero frequency and large k, because the Markov approximation becomes exact at C = o. Noting that C44(k)/p o (l/Pm) (cf. Equation (3.5)) we obtain Rj(k,o) from (3.13) as (2pm) /k whereas 1/2 (3.17) yields (Hpm/2) / /k. The ratio is 74/H - 1.12. Thus (3.13) recaptures the ideal gas result in the limit of small w and large k. It may be pointed out here that (3.10) can be approximated for large frequencies by replacing cp22(k,i1) by P22(k,o)/ic (short-time expansion of 22(k,t)). Since our inter13

est lies in the small frequency region we shall not dwell on this point further even though c22(k,o) is calculable exactly. With the aid of (5.20) we have been able to obtain an expression for the transverse current power spectral density (3513), which contains only one adjustable parameter, k. The value of k determines the transition from the o o small to large k limits. It is expected to be in the vicinity of the main peak of the structure function S(k), which occurs at k = 2A-1 in Rahman's computer (2) calculations for liquid argon. We have chosen k = 1.5'-1 which yielded 0 the best fit to the computed curves although the value of k is found to be not too critical. The other constant in (3520) is w2(o) which is obtained from (5.15) as Ws(o) =.1 x 1026 sec, corresponding to a value for the shear viscosity o = 2.8 x 10-3 poise at p = 1.407 g/cm3 and T = 76~K for liquid argon. The values of C44(k) were computed according to Eq. (3.5) in which the interparticle potential is taken, following Rahman, as V(R) 4 6=( with (E/kB) = 120~K (kB = Boltzmann constant) and a = 3.9 A. The variation of C44(k) is plotted in Fig. 1. Figures 2 and 3 show the variation of sg(k) and Ts(k,o) = C44(k)/0S(k) with k. We observe that the k-dependent shear viscosity decreases very rapidly by a factor of 100 in the region of k from zero to 2 A-1, and approaches zero as (l/k). This k-dependence of r.s(k,o) appears to be crucial to the behavior of Rl(kCw) for the k-values in 1-4 A1'. Figures 4 and 5 show comparison between the calculated. curves and Rahman's 14

(2) data. ) It is noteworthy that the present model predicts well the cut-off wavelength in the dispersion curve, i.e., wsm(k). Other features are self-explanatory. 15

4. A GENERALIZED HYDRODYNAMIC DESCRIPTION The previous application indicates that the correlation function associated with a dynamical variable al(t) (e.g., Jl(k,t)) can be obtained by solving the appropriate generalized. Langevin equation. If only the autocorrelation function is of interest, the one-component description of the system is sufficient in principle~ The correlation function in this case is obtained by solving Rl(t) + f du cp1(t-u)Ri(u) = 0 (4.1) (Note that the frequency matrix is always zero in one-dimensional description.) The damping function pl(t) involves f(t) = exp[t(l-P)iL](l-P)al(o) where P projects a phase function onto al(o). Although (4.1) is exact, the calculation of cpl(t) is as difficult as calculating <al(t)al> directly. Crude approximations for cpl, such as the Markov assumption, are generally not precise enough to include even the qualitative features of the correlation function, or the power spectral density associated with it, for large values of w and k. By introducing instead, a multidimensional description of the system, one actually extracts a great deal of information about the collective motion of the system through the frequency matrix even though one may still be interested only in the autocorrelation function of a single variable. This information is contained in pl(t) in one-dimensional description. A proper choice of the additional variables in a given system, can lead to a sufficiently precise expression for the correlation function in a wide range of w,k even with crude 16

approximations on the multidimensional damping function. The variables Jl(k) and H113(k) introduced in Section 3 provide such a description for the transverse-current correlations. The description of the longitudinal current correlations requires a more detailed description of the fluid including thermal and viscosity effects as will be demonstrated in this section. Moreover, a multidimensional description also allows the computation of the various crosscorrelations between the variables in the set in terms of the same thermodynamic and transport parameters. The purpose of this section is then to present a. 14-dimensional description of a simple liquid, and to compute specifically the transverse and longitudinal current correlations. This description includes thermal effects, and sheds light on the anisotropies in the fluid. for large k values. For state variables, choose a col[p,O,a,Ji, q], (4.2a) where J. and q. are vectors with three components (i,j = 1,2,3) and a is a 6component vector with. = 1,..6, viz., a = col[o11, 022, 33, a13 G23, 13]. (4.2b) The variables 0, a and q. are defined by N k ikox p(k) = m E e (4.3) a=l 9(k) = E(k) - )p- p(k] (4.4) <p(k) p (k)> 17

<1.. )p (k)k)> < (k)(k)> o. (k) -I..(k) - -p(k) - - - (k) (4.5) j 1J <p(k)p (k)> <0(k)9 (k)> q.(k) Q.(k) - <Q J (k)>.<J(k)J*(k)>-'lJ(k) (4.6) The definitions of J(k) and HI..(k) have already been given in (3.2a) and (3.2b) 1o - respectively. The quantities E(k) and Q.(k) are the energy density, and:the energy current density respectively. They are defined as N N E(k) 1 a + 1 mv V(I x_ a P| ik-xx e4. x Q.(k) = 1 [ my va+ 1 L V(xC-x eik x 223: lX 2 The following usual conservation laws prevail among the variables p, J., H1.. E, and Q.:'.ap1\6 -+ ikQ. P ) (4.9c) 18ik, (4.9b) t - 18

The tensors HI. and.. are symmetric and have only six independent components as implied in (4.2b). The variables aij and qj denote the viscous stress tensor and the heat flux vector. The 9(t), defined by (4.4), will be replaced later by 0(k) T(k) (4.10) p C (k) o v whose average with respect to a perturbed distribution function yields the tem(l6) perature in the conventional linearized hydrodynamic description of a fluid. (16) Such an identification, however, is not needed for the presento The quantity C (k) will be defined later (it will be identified as the specific heat at constant volume). The average values of p(k), E(k), and HI. (k) are zero for k + 0. When (21) k = 0, we have' <p(0)> = p (4.11a) V 0 P P 2 <E(0)> -= - + - I d3R)(V(R) (g(R)) (411b) V 2 m 2 m2 - 2 1 Io i ^ t0^ ~ ^j ^ = -Po P o 0 SR dV(R) V ij- (0)i o d g 6) d R g(R (4.11c) Here, P is the equilibrium pressure. We shall always assume that these averages are subtracted from p, E, and HII. whenever they are not zero, so that the state vector a will denote the fluctuations. Eight of the 14 components of the state vector a are even, and the remain19

ing six are odd functions of particle velocities. Hence a can be decomposed into even and odd parts as a = col[a,0] + col[O,a]. Consequently, the static correlation matrix _ = <aa > splits into two disjoint submatrices as 4e 0 r = ~ (4.12) where e e e* *<PP > 0 0 <22 1 and 0 1 <JJ*> 0 ~ =<a a~ > = q | q (4.13b) L 0 <qq > e O The block-diagonality of 0 and. is a consequence of the choice of the state variables as in (4.4), (4.5), and (4.6), which imply the following orthogonality relations: <p > = <a..p > = <. e > = <q.J> = 0. (4.14) The static correlation functions appearing in (4.13a) and (4.13b) will be discussed latero The frequency matrix Q = <a a >of- can be written as 20

.e o* o01 0 <<a a >. iQ = ~1 _ _ _ |(4.15) e*o e- <a, >_ 0 As a consequence of the conservation relations (4.9), f(o) = (l-P)a has only nine nonzero components: f(o) = col[0,O, fO,O ]. (4.16) Therefore, the damping matrix Tp(t) is of the following form: 0 0 0 0 0o ar(t) 0 a(t) iT(t) =, _ (4.17) a0 0 0 Lo q-(t) 0 q(t where p —and p -- are 6 x 6 and 3 x 3 square matrices. The off-diagonal matrices cp —(t) and c_- are 3 x 6 and, 6 x 3. Substituting (4.15) and (4.17) into the generalized Langevin equation discussed in Section 2, we obtain the following set of equations: p(t) = ik.J(t) (4.18) <Jp > <Je > J(t) -(t - e(t) = ik.a(t) (4.19) <pp > <Ge > e(t) + <J > ><JJ >-l J(t) = ik-q(t) (4.20) 21

b0 * Jt -1 t Crq (t) <q ><qq > q(t) + / _-(t-u)'q(u)du - <*J <JJ >.J(t) + S — (t-u) a(u)du = f(t), (4.21) q(4t) - - e(t) + f P —(t-u) q(u) <00" >'_* * -1 qt q q - <q >'<a > oe_(t) + f -— (t-u)'a(u)du = f(t) (4.22) n writing (4.20) and (4.22) we have used the orthogonality properties of the variables, and the identity: <qp > = 0 (4.23) These follow from the sum rule <AB > = -<AB > where A and B are two arbitrary (21) dynamical variables. ) Notice that we recover the conservation laws (4.9) in Eqs. (4.18), (4.19), and (4.20). Equations (4.18) through (4.22) describe the time evolution of the state vector a = col[p,e,a,Ji,q.] exactly. We shall approximate them by neglecting the coupling between the viscous stress tensor and heat flux vector and introducing the Markov assumption in (4.21) and (4.22). The first approximation is equivalent to setting <aq > 0 and -— (t) 0 in (4.21) and (4.22), and the second to t —(t-u).a(u)du W-(k).a(t) (4.24a) _q-(t-u)q(u)du Wq(k).q(t) (4.24b) 22

where wr 0f a Cr W (k) S dt — (t), (4.25a) w0(k) = C dt q(t) (4.25b) Furthermore, the following definitions and equalities will be introduced: * <Jp > CL(k)ik (4.26a) Lt L <PP > <JT > - * CL(k)p L(k)ik (4.26b) <TT > 1 * P - <i.J.> -- j (4.26c) V i j P ij -<ee >.>Ep <0 0>.= V L<EE > I T C (k) (4.26d) <JiJl > oL <pp > * * <qe > * i.k KT-(k) (4.26e) <e > < a > - > (k) (4.26f) V pt v 4jv Then, we obtain the following approximate description: p(t) = ik.J.(t) (4.27a) j.(t) - Ca(k)ikj[p(t)+L(k)(pT(t)] = ik am (t) (4.27b) 23

0(t) - 5L(k)C (k)T p(t) = ik.qj(t) (4.27c) LL i o j (t) - Z (k)e (t) + W (k)o (t) = f-(t) (4.27d) qj(t) - K. (k)ikj (t) + W q (t) fJ(t) (4.27e) L m t jmm J where ~ is the rate of strain tensor, i.e., mn (k) = [k J +k J ] (4.28) mn 2p mn nm o We shall discuss the physical implications of the various quantities appearing in this set. The C2(k), introduced in (4.26a) can be defined also by L C2(k) =. (4.29a) L * <pp >' i~mS (k) ~(4.29b) funS(k) where S(k) is the structure factor defined by S(k) = 1 + / d3Re i- [g(R)-il], (4.29c) m g(R) being the static pair correlation function. We will refer to C (k) as L the longitudinal isothermal speed of sound. Its relation to the conventional isothermal speed C (k) can be established by considering the definition of the (16)(21) latter: <(Tr i )P >. 13 = 24(k) (4?29d) <pp > T 24

(16,21) where P is the thermodynamic part of pressure) p = (1/3)Trll. (the nonij thermodynamic part, by definition, has no projection on p and T). If we define a transverse isothermal speed as <22p > _ <1p > - c2(k) (4.o0) *- * 2L <PP > <PP > then, C2= (1/3)(CL+Cq). In the isotropic limit where ka << 1, a being the o L mean linear force range, then 2 C = = /2)C2, and the distinction between o L L longitudinal and transverse isothermal speeds becomes irrelevant. It is important to note that only longitudinal isothermal speed. is related. to the structure function S(k) as in (4.29b). The quantity 5L(k) in (4o26b), which we refer to as the longitudinal thermal expansion coefficient, may be equivalently defined by < * - C2(k)p Po(k). (43.1) *L oL <TT > We may also define a transverse expansion coefficient pT(k), similar to (4.30), as (p oTC /2) = <1T >/<TT > = <122T >/<TT >. The conventional definition of P(k) as a derivative of the thermodynamic part of pressure is <Tri. T > P C2(k)po3 (k) * aT -(4.32) <TT > Then, p is related to pL and pT by LC L2+TCT o = -- * (4.33) L T 25

In the isotropic limit, P = PL = PT. Equation (4.26d) defines the specific heat at constant volume. This de(21) finition is identical to that by Schofield.( The quantity KL introduced in (5.26e) can be equivalently defined by * - 1K11*i (<.5>) KL(k) - <q > >1 (4 34). =L* 0 *QsQs> - V' <ee > <9O > o It will be related later to the components of the thermal diffusivity tensor (cf. 4.42). The symmetric tensor Z (k) defined by (4.26f) is related to the elastic 4V moduli tensor C by fIv <n p X>pH > <1 T XTH >I E c= - c + * (4.35a) V 4 <PP > <TT > where C <[ nI > (43.5b) 4V V 4 v We have already encountered some elements of C. In (355) we have used C44 = (P/V)<1T31I31>. Since <T3jp > = Fj3<s33P > and <l3jT > = 5j3<n33T >, the terms in parenthesis are zero. We then obtain'44 = C44 = C55 = 5 (4.36) we shall also need 733 and Z32, in terms of C33 -= /V <T3I33> and C32 = P/V <n3lT2>(clearly E32 = Z31, and C32 = C31). From (4.35a) we get 733 = C33 - po y(k)C(k) (4.37a) 26

where we have introduced P2(k)C2(k) v(k) r + C2k) T]. (4.37b) (a consistent notation would be Y33. For simplicity we let y = 733). The quantity 7(k) reduces to the conventional ratio y = (C /C ) in the limit of k - 0 (the product 47 CL(k) can be interpreted as the longitudinal adiabatic speed. of sound. However, such identifications for large values of k become ambiguous, and in fact are not needed. The important point is that all the elements of Z are expressible in terms of the integrals involving the two[IV and three-particle static correlation function and interpartical potential.) For 732 w'e find from (4h35a) 32 = C32 - P C(k)732 (4.38a) o L where <piiHg> > <Tn > /, p 732 = + PL < (4.38b) <J13i> <JlJil> We note that 732 is different from y in (4-37b). They become identical in the isotropic limit mentioned above. It is known that in a crystal with cubical symmetry, the elastic moduli (23) tensor has only three independent components as discussed by Kittel: (23) ClI C12 C12 C12 Cll C12 0 c C12 C12 C11 (4.39) ][ IV0 C44 0 (4 C44 0 0 C44 27

As an approximation we shall assume that both Z and C have this form. The ~v ~v numerical values of C11, C12 and C44 are plotted in Fig. 1. The values of CL(k) L are determined:from the knowledge of S(k) (cf. 4.29b). Therefore, if we approximate y and 732 by their limit for small k, we can compute the elements Z12, 11 and Z44. We now focus our attention upon (4.27d) and (4.27e) which are generalizations of Maxwell's model for relaxation of the viscous stress tensor and Fourier's law. They both contain the relaxation frequency matrices W- and 1-v Wq. Using the former and the relation a. = E valid for small frejm 1ij ijm Im quencies wave numbers, we may define a "viscosity" tensor T as z = W- (4.4o) [v rI4 Wv It is consistent with (4.39) to assume the same symmetry for n. as Z. Then, the viscosity tensor also contains only three independent elements q11, 112 and 144o Since the inverse of a matrix of the form of (4.39) has also the same form we find that the frequency matrix W posses only three independent frequencies c11, CD12 and C44: XCo = -11111j2(Z1j-2-12) (4.41a) (Tll-12)(T1-+2T12) W12 - (111-112)(\1111212) (4.41b) C.44 = (4.41c) 144 These frequencies as well as components of the viscosity tensor are functions of k, which have to be determined, from either (4o25a), or from the asymptotic 28

behavior of the correlation functions as discussed below. In the case of Wq., we define a thermal diffusivity DT(k) D (k) = [W]. KL(k) (4.42a) On the basis of symmetry arguments it is reasonable to assume that W-. in q q q la diagonal with elements ll1 = 22 and W33 - o. Only X appears in the calculation of the longitudinal current correlations. The corresponding component of D is _ (k) - KL(k) DT(k w ) (4.42b) 0 which reduces to the conventional thermal diffusivity in the limit of k -* 0. TRANSVERSE EQUATIONS The projection of J on a fixed axis perpendicular to k satisfies the following equations obtained, from (4.19) and (4.20) Jl(t) = ikG4(t), (4.43a) ik ( 4(t) - C44Jl(t) + W ( = (t) = f(t) (4.43b) P 4 4V V The assumption of cubical symmetries removes the coupling between o4(t) and other components of a (t), because it implies W- = 5. If we let v 4v 4v 44 CUl = Ws in (4.43), they become identical to (3.8) in the Markov limit, which was obtained by using a two-component description. Hence, we can use (3.14) to compute L44(k)o Using (30llb) and (4.41c) we also find that 044(k) = Ts(k,O) (4.44) 29

which enables us to interpret p44(k) as the k-dependent shear viscosity. It is important to note that Eq. (3.8), corresponding to (4.43) in the hydrodynamic description, is exact; whereas (4.43) are approximate. It therefore appears that the damping term in (5.8b) automatically includes all the coupling effects which are ignored in (4.43b). The definition of )44 obtained from (4.25a) as a component of the multidimensional damping matrix is not identical to (3.14) because they involve different projection operators. LOGITUDINAL EQUATIONS The compenent of J parallel to k satisfies the following set of equations: p(t) = ikJ3(t) J3(t) - ikCL[p(t)+LoT(t)] = ika3(t) 3(t) - Z11J3(t) W11(t) +2(t+ [1103(t) = (t) P 3 al(t) - - Z12J3(t) + W12C(t) + W11C2(t) + )12a3(t) = (t) (4.45) p 1 2(t) - i E12J3(t) + W11o(t) + c1202(t) + W12a3(t) (t) 0 2 @(t) - LC2To(t) = ikq3(t) 43(t) +- o[q3(t)-ikD 9(t)] = fa(t) T 3 -j where we have already introduced the cubical symmetry. Multiplying these equations by J3(o0) taking ensemble averages, using <f-(t)J3(o)> = <f(t)J3(o)> = 0, and Laplace transforming with respect to time, we obtain the following with rese following expression 50

for the longitudinal current power spectral density Rl (k,o) = Jdt cosact Js(t)J(46a) 0 * (4.46a) <J3J3> as /lk u L k kC 2 R(k) = i - L+ (kiW) + -)CLk (7-1) 1 (4.46b) (kFc L o i na+k2 ( k. i )j where p * L (k, ico).<3 (4.47a) <J3 (icu)j3> iR(k,i ) + 5 rl'(k icu) (4.47b) 3 Po <1/3 Troi. (i) J3> + B (k, ic) i- *+ (4.48a) B'ik <J3(iW)J3> 5 JiJ+( B+ 22 1+21 (4.48b) "B = i21+2%d2 =48) s(kic) 2 ic+c' 4.4a C - CL1 - C2 = (4.49b) 6(k, i) <q(i)J> (450a) ik<o( i ) J3> DT (k) l +i ((l/Ot)^O) (4. 50b) 31

The variations of Rll(k,w) with k and C has also been computed by Rahman(2) for an argon-like simple liquid. We shall interpret his data using an approximate form of (4.46). It is important to note that (4.46) is an exact result if TL(kiw) and 5(k,iuc) can be evaluated exactly using their definitions in (4.47a) and (4,50a) in terms of correlation function. The approximations which we have introduced above (i.e., the neglect of the coupling between heat current and viscous stress tensor, the Markov assumption, and the assumption of cubical symmetry), were actually made to evaluate L (k,iw) and 5(k,ic) explicitly. We shall discuss them in some detail. We observe in (4o47a) that TL(kin) is the sum of two terms. The first term, B (k,iw) may be identified as a generalization of the bulk viscosity as its definition in (4.48a) indicates. The explicit form in (4.48) has been obtained from (4.45) considering Troij = aC+a2+a. The second term Ts(k,iw) defined by (4.49a) is not the k- and a-dependent shear viscosity qs(kic) (cf. 3o11) as one might have expected. For small k limit, however, this is indeed the case. In the limit of k - 0, the components of (21,23) the elastic moduli tensor satisfies ( ) C11 - C12 = 2C44. (4.51) Since y = 732 in this limit, we also have Z1 —Z12 = 2Z44 = 2C44. Furthermore, %1-T1 = 42 244 because n and Z are assumed to have the same form. Then, pIv [iv (4.49) yields cll1-k2 = (44, ioe., the relaxation frequencies for the viscous stress tensor are not independent. Substituting these in (4.49) we obtain sI(k,iw) -+ s(k,iw) (isotropic limit) (4.52a) and 32

L(k,io) B+ q(kiu) + (4) S(k iwu) (4.52b) If, as an additional approximation we assume Uo >> W(12 then, aB w C =c S 111 CL and we have for all k Cll(k)-p 7(k)C2(k) \L(k' i).it() (4. 53) Physically, this assumption is equivalent to ignoring the coupling between 03 and oa and ~2 in (4.45) (set ac2 = 0). An other consequence of this approximation is that the transverse relaxation frequency us = 044 approaches cL for small k because d44 -+ Cj11-C2 L C01. These observations have been verified approximately using Rahman's computer results as will shortly be discussed. As a result of the k2 factor multiplying 6(k,icW) in (4.46b), the thermal effects become insignificant for large values of k (e.g., k > 1A 1 for argon), and Rll(k, e) reduces to. C1L k2 Rll(k,w) = Re -e- + - pL(ki2) (4.54) We have used (4.4) and (43) to interpret Rahman's data (we have set y(k) 1 We have used (4-.54) and (4.55) to interpret Rahman's data (we have set 7(k) - 1 as a further numerical simplification) which may be called the isothermal approximation: 2 (k)k2[Cll(k)-po CL(k)]/po Rll(k, w) L L. (4.) 2(k) [W2-k2C2(k) ] 2+(2-kCkC (k)/p)2 L L The only unknown in this expression is nL(k). Its value for k - 0 was 33

estimated using CLt)-PC2(o) Cll(~)-PoCL( (456) WL(0) = (4. 6) iB+3 s where TB and s are the conventional bulk and shear viscosities, respectively. We note that lim C11(k) = K + G (4.57) 00 3 00 k->oo where G was defined in Section 3, and K is the bulk modulus in Zwanzig's nota00 00 tion.( The value of C1(o) was read from Fig. 1, and that of C2(o) (the 1L ordinary isothermal speed of sound) was calculated from (4.29b) (Rahman has also calculated S(k)). To model the large k behavior of w (k) we consider the ideal gas, L IG2 2 Rll(k, ) = Pm - - e- pm2/2k (4.58) k2 2k2 which has a maximum at C2(k) = 2k2/pm. Requiring that the frequency at which (4.55) m attains a maximumfor a fixed k approach2k2/Pm for large k, and using the same interpolation formula as in (3.20) we obtain the following expression for LL(k): L(k) =8k (Cll(k)-po CL(k)-Po/pm)/p oL(k) = 82 2 c2(o) k2(C)-(k) k) )-2po/m)/p +., (4.59) [l+k2/k2] where w L(O) is given by (4.56). The variation of wL(k) is shown inFig. 2 34

(k = 1.5 A as in the case of transverse correlations). IG It is interesting to compare again Rll(k, )/2a from (4.13a) and Rll(k, )/o from (4.18) for k -+ co and +- 0, as in the case of transverse current power spectral density. Using C2(k) - (l/pm) and Cll(k) + (3P /pm) as k -+, we find -2 3/2 r-IG Rll(ko)/o + (Pm)3/2 /2 /k3 and R11(k,o)/f 2 (m)3 /2 X/i /k3. Thus, the approximate formula for Rl(k,cb) yields the ideal gas limit for small X and large k within a factor \T3/IT 0.98. The curves in Figso 7 and 8 are calculated using (4o55) and (4o59)o The agreement between Rahman's computer results and the theoretical curves are remarkably goodo The thermal effects dominate for small values of k such as those involved in light scattering. We then approximate (4.46b) as / k22(o 2 C (o)k2 (y-) CL k2 Rll(k, c) = Rei( - +) +. (4.6o) C a L iw~+k2D (o ) -which is the result corresponding to the conventional hydrodynamic description with constant transport parameters. Since the variation of the thermal diffusivity DT(k), thermal relaxation frequency and y(k) are not readily available, we may use their limiting values as k - 0 in (4.46b), and substitute TL(k,co) from (4o59) to find R k(k, C) = k Cll (k)-Po CL(k) Rll(k, ) = Reid -() - -.O L y l+(0 ) 35) kD (o) 1 +kC(k)Q(-l) i+ T.j (4.61) 55

where w (k) is to be taken from (4.59). All the functions of k appearing in L this expression are calculable in terms of S(k) (or g(R)) and interparticle potential, and the constants 7, DT(o) qL(oo) are experimentally available. CONCLUSIONS In this paper we have used the generalized Langevin equation to investigate fluctuations in simple liquids in the frame-work of a generalized hydrodynamic descriptiono The various transport coefficients appearing in this description have been defined in terms of time-correlations of the dynamical variables, and computed numerically whenever possible for argon-like liquid.s The results, such as the variation of viscosities with wavelength, may be applicable to other simple liquids. In view of the good quantitative agreement with Rahman's data for all values of k and CX encountered in neutron scattering, we may expect the final expression (4o61) to be applicable to the interpretation of coherent neutron scattering as well as light scattering from dense fluids. (Figure 8 represents W2S(k,w) = R1(k,w) as a function of k and w. No attempts have been made to compare these to actual neutron scattering data, because the discussion of the contribution of incoherent scattering which is important in natural argon'( is beyond the scope of this paper.) The formalism developed in this paper enables one to calculate correlations between other pairs of hydrodynamic variables, eog., <E(t)E > and <Q(t)Q >, <Jl(t)J11>, etc. Computer results for these correlation functions in conjunction with the analytical calculations will shed light on the k- and cu-dependence of the thermal parameters y(k), W (k) and DT(k) as we have demonstrated for the viscosities. 36

ACKNOWLEDGMENTS We wish to thank Professors R. K. Osborn and. J. Duderstadt for their helpful comments and reviewing and editing of this paper. Also we thank Professor S. Yip for discussions and providing his own related calculations prior to publication. The work described above was supported. in part by the National Science Foundation. 37

REFERENCES (1) Rahman, A., Phys. Rev., 136, A405 (1964). Rahman, A., Phys. Rev. Letters, 19, 420 (1967). (2) Rahman, A., in Neutron Inelastic Scattering (International Atomic Energy Agency, Vienna, 1968), Vol. I. (3) Nelkin, M, VanLeeuwen, J.M.Jo, and Yip, S., Inelastic Scattering of Neutrons in Solids and Liquids, Vienna (1965), Vol. II, 35-57(4) VanLeeuwen, J.M.J. and Yip, S., Phys. Rev. 139, A1138 (1965). (5) Nelkin, M. and Yip, S., Phys. Fluids, 9, 380 (1966). (6) Greytak, Jo J, and Benedek, G. B., Phys. Rev., Letters, 17, 179 (1966). (7) Nelkin, M. and Ranganathan, S., Phys. Rev., 164, 222 (1967), Ranganathan, S. and Nelkin, M., Journ. Chem. Phys., 47, 4056 (1967). (8) Nelkin, M. and Ortoleva, P. J., IAEA Symp. Inelastic Neutron Scattering (Copenhagen, 1969), (9) Callen H. B. and Welton, T. A., Phys. Rev., 83, 34 (1951), Kubo, R., Tokyo Summer Lectures in Theoretical Physics, Part 1, Many-Body Theory (1966). (10) Landaue, L. D. and Litschitz, E. M., Statistical Physics (Pergamon Press, London-Paris, 1959). (11) Rytov, S. M., Soviet Phys. JETP, 6, 130 (1958). (12) R. Mountain, Rev. Mod. Phys., 38, 205 (1966). (13) Kadonoff Lo P. and Martin, P. C., Annals of Physics, 24, 419 (1965). (14) Chung, C. H. and Yip, So, to be published in Phys. Revo (1969), 38

REFERENCES (Concluded) (15) Zwanzig, R., "Statistical mechanics of irreversibility," in Lectures in Theoretical Physics, edo by W. E. Brittin, (John Wiley, New York, 1961), Vol. 3. (16) Mori, H., Progress of Theoretical Physics, Vol. 33, No. 3, March 1965. (17) Mori, H., Progress of Theoretical Physics, 34, 399 (1965). (18) Akcasu, A. Z. and Duderstadt, J. J., to be published in Phys. Rev. (19) Zwanzig, R., Phys. RevO, 144, 170 (1966). (20) Wang, M. C. and Uhlenbeck, G. E., Rev. Mod. Phys., 17, 323 (1965); Chandrasekhar, S., Rev. Mod, Phys., 15, 1 (1943). (21) Schofield, P., Proco Physo Soco, 88, 149 (1966). (22) Zwanzig, R., and Mountain, R., J. Chem. Phys., 43, 4464 (1965). (23) Kittel, Co, Introduction to Solid State Physics (John Wiley, New York, 1953). 39

FIGURE CAPTIONS Fig. 1. Variation of elastic moduli with wavelength. Fig. 2. The variation of the transverse and longitudinal relaxation frequencies with wavelengths. Fig. 35 The variation of shear viscosity with wavelength. Fig. 4o The transverse current-current correlation, Rl(k,o), vso. c for various k. Fig. 5.1. The maxima of a transverse current-current correlation function vs. wave number, k. Solid line due to present theory and points represent Rahman's data. Fig. 5.2. The frequency at which the transverse current-current correlation function is maximum. Fig. 5~35 The transverse current-current correlation at zero frequency. Fig. 6. The longitudinal current-current correlation function vs. C for various ko The points are from Rahman. Figo 7.1. The maxima of the longitudinal current-current correlation function vs. wave number, k. The points are from Rahman. Fig. 7.2o Frequency at which the longitudinal current-current correlation function is maximum. Fig. 8. The scattering function, S(k,w), vs. c for various k. 40

I00 80 60 40- \* —C*03/6)Cll 80t ^ ^4/3 C44-^ t(4/31'3/ 2- ((-3/E)Cl2 C12 " 20 8I ~2 k iA 5 6 7 k in. AFig. 1 41

1.6 ~1.4 WL (k): LONGITUDINAL 1.2 1.0 U _. - Os (k): TRANSVERSE a.6.2 0 1.0 2.0 3.0 4.0 5.0 I-I k in A Fig. 2 42

3 2 1.0.8 w.6 (c 0 -..4 -J -- \ z 0.2.08.06.04-.02 0 1.0 2.0 3.0 4.0 5.0 k in A' Fig. 3 43

2.4 2.2 2.0 K.3-9 / n CALCULATED 1.6'o 1.4 - z 1.2 K=3.7AC00MPUTER 1.0 =- 0.8 0.6 0.4 0.23.6 3.44 3.2 3.0 3.8 K= 0.5 A-I W 2.6 2.4 z 2.2 2.0 1.8 1.6 K=2.5 A0 1.4 K=.3A1.2 1.0 0.8 0.6 -K 3.7 0.4 0.2 0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0 0 0.2 0.4 0.6 08 1.0 1.2 1.4I 1.6 1.8 20 wIN 10'S SECI wIN 10'S SEC.I Fig. 4 44

5.0- Fig. 5.1 U 4.0-'0 c 303 \ a 2.0 o 1.0.-.....~..~ ee 0 i I I, l l l 0.4.8 1.2 1.6 20 2.4 2.8 32 K in -I Fig. 5.2.5 o.43".2.1 0.4.8 1.2 1.6 2.0 2.4 2S 3.2 K inA5.0- \Fig. 5.3 - 4.03.0 - c 2.0 ~ *O 0 00.0.000 0 1.0 O 4.8 1.2 1.6 2.0 2.4 2.8 3.2 K in A45

3.6 - 3.6 - 3.6 - 3.2 3.2 - 3.22.8 2.8 - 2.8 - 24 K= -5$' 2.4 - 2.4<A K=.9$-' K= 1.7?-I 2.0 2.0/- 2 01.6 - 1.6 - 1.6 1.2 -.2 - 1.2 - 2.8 -.8 -.8.4.4.4 20.C - 2.0 - K 0-^S = 0.4.8 1.2 1.6 2.0 0.4.8 1.2 1.6 2.0 0.4.8 1.2 1.6 2.0 3.6 3.6 - 3.6 3.2- 3.2 - 1 3.2 2.8 2.8 - 2.8 2.4 0 2 2.4K= 2.1 24' 1.6 - 1.6.8 In 1.60. 1.2 - 1.2 - 1.2.4.8 1.2 1.6 2.0 0.4.8 1.2 1.6 2.0 0.4.8 1.2 1.6 2.0 cd in 10'3 sec. Fig. 6 46

5.0 0 4.0c _ 3.0* Ile ~~~ 1.0Fig, 7.1 o I I I i 0 4.8 1.2 1.6 2.0 24 28 3.2 3.6 K in -I.8t- ~.4.8 1. 2 1.6 2.0 2.4 2.8.3.2 3.6 47 4~7

100lo 10 _ 1.0 80 8.8 60 6-.6 40\ 4.4 K = 2.9' 206 2 -.2 - 10 1~~~~~~~~~~~~~~~~~~~~.0 I.4 -.0 - \.00 -.8 K23.3 0821 6*.6.06 34 K=2.1A.4.04 2-.2 -02 K 2.5'A 1.0 K=1.7A.c i.01.8 -.08 -.008.6.06 -.006.4.04 -.004.2 -.02 -.002.1,.01, -.001, 0.2.4.6.8 1.0 1.2 0.2.4.6.8 1.0 1.2 0.2.4 6 8.0 1.2 w in 10'1 sec-1 w in 10'3sec-I w in 13 sec1 Fig. 8