NONEQUILIBRIUM FLUCTUATIONS AND TRANSPORT
IN SHEARED FLUIDS
BY
JAMES FRANCIS LUTSKO
A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL
OF THE UNIVERSITY OF FLORIDA
IN PARTIAL FULFILLMENT OF THE REQUIREMENTS
FOR THE DEGREE OF DOCTOR OF PHILOSOPHY
UNIVERSITY OF FLORIDA 1986
To my parents.
ACKNOWLEDGMENTS
I would like to thank Jim Dufty for his guidance, encouragement and patience throughout the time I have worked with him. His supervision and friendship have proven invaluable in the course of my work.
In addition, my gratitude goes to the other faculty who have guided me through my studies as well as to my fellow travelers in graduate school. In random order, they are Dave, Brad, Pradeep, Gary, Jim, Simon, Tom, Bob and Marti. Thanks also go to Bob Caldwell for assistance with the numerical analysis. Finally, special thanks are extended to Stephanie for her support, encouragement and friendship.
iii
TABLE OF CONTENTS
Page
ACKNOWLEDGMENTS..................................................iii
ABSTRACT........................................................... vi
CHAPTERS
1 INTRODUCTION................................................1
2 PHENOMENOLOGY AND SIMULATION .............................12
Phenomenology...........................................12
Experiments ...............................................13
Simulations ...............................................15
3 FLUCTUATION MODELS......................................22
Generalized Langevin Equation ............................22
Long Wavelength Fluctuations .............................29
Small Length Scale Fluctuations ..........................32
4 LINEARIZED HYDRODYNAMIC MODEL.............................. 34
Linearized Model........................................... 34
Green's Function............................ .............37
Correlation Functions.................................... 47
Physical Properties......................................52
Probability Densities................ .................57
5 NONLINEAR FLUCTUATING HYDRODYNAMICS........................63
Naive Renormalization ........................................63
Formal Renormalization........................... ........69
Application of Formal Renormalization ..................... 74
6 SHORT TIME MODEL........................................80
General Properties of the Reduced Distribution
Functions ...............................................80
Explicit Form of Matrix Elements..........................83
iv
Approximations ............. .............................86
Hydrodynamic Equations..................................90
7 ANALYSIS OF THE SHORT TIME MODEL..........................96
Equilibrium Generalized Hydrodynamics.....................96
Nonequilibrium Generalized Hydrodynamics..................97
8 CONCLUSIONS........................................... 121
APPENDICES
A LINEARIZATION OF THE LANGEVIN EQUATION ..................125
B GENERALIZED EIGENVALUE PROBLEM ..........................133
C DYNAMIC STRUCTURE FACTOR ................................139
D ELIMINATION OF THE THREE POINT FUNCTION .................142
BIBLIOGRAPHY......... ....... ......................................46
BIOGRAPHICAL SKETCH..............................................150
V
Abstract of Dissertation Presented to the Graduate School of the University of Florida in Partial Fulfillment of the
Requirements for the Degree of Doctor of Philosophy
NONEQUILIBRIUM FLUCTUATIONS AND TRANSPORT IN SHEARED FLUIDS
BY
JAMES FRANCIS LUTSKO
December 1986
Chairman: James W. Dufty
Major Department: Physics
Hydrodynamic fluctuations in a sheared fluid far from equilibrium are described by means of nonlinear stochastic differential equations. Two models resulting from two different limits of the exact FokkerPlanck equation describing hydrodynamic fluctuations in an arbitrary fluid are presented. The first model results from restricting attention to fluctuations which take place on length scales large with respect to the mean free path. The resulting NavierStokesLangevin equations are specialized to shear flow and the associated linear hydrodynamic modes are determined. These are used to calculate the correlation functions of fluctuations around the steady state with the full nonlinear dependence on shear rate necessary for large shear rates. It is then shown how the coupling of the linear modes through the nonlinearities in the stochastic equations allows one to calculate the effect of the nonequilibrium state on the macroscopic transport properties.
vi
The second model results from the short time limit of the exact FokkerPlanck equation and describes fluctuations even on atomic length scales provided the mean free path is small enough. The linear modes of this model are obtained and it is found that they differ considerably from the long wavelength modes. In particular, it is found that at large shear there are two sets of propagating or sound modes. A large shear instability is also found in this model and is tentatively identified with a large shear dynamic phase transition observed in computer simulations of sheared fluids.
vii
CHAPTER 1
INTRODUCTION
Statistical mechanics is a methodology for determining the
physical properties of many body systems from the microscopic dynamics of its constituents. It is applicable when the properties one wishes to determine arise as a result of the interaction of a large number of the elementary constituents or when the time scales of interest are large with respect to the typical interaction time of the constituents. When this is the case, attention is shifted from the dynamics of individual constituents to the density in phase space (i.e., the space defined by the degrees of freedom of the system) of the system as a whole which is characterized by the distribution function of the system. The dynamics of the distribution function is given by the Liouville equation, which expresses the conservation of phase space density in a conservative system, together with appropriate boundary conditions, initial conditions and source terms representing the action of nonconservative external interactions.
When the distribution function is constant in time, the system is said to be in a steady state. Equilibrium is the particular steady state obtained by conservative systems and, in general, is unique in that the associated distribution function is known. For this reason, the statistical mechanics of equilibrium systems is well developed and
1
2
quite successful at predicting macroscopic properties such as transport equations and thermodynamic susceptibilities.
The statistical mechanics of systems not in equilibrium is, in contrast, poorly understood. This difficulty arises from the fact that nonequilibrium states are difficult to characterize; the distribution function is not known. For systems close to equilibrium linear response theory provides a means of studying nonequilibrium systems by relating their properties to those of equilibrium systems (e.g., fluctuationdissipation theorems relate transport properties to equilibrium fluctuations). However, far from equilibrium one has only the Liouville equation.
To get a handle on systems far from equilibrium, one typically restricts attention to some subset of the degrees of freedom which obeys an approximately closed dynamics on the length and time scales of interest. Closed, in this context,.means that couplings to the degrees of freedom not in the subset being studied are negligible on the average. For example, the macroscopic equations of continuum mechanics for the local conserved densities (e.g., hydrodynamics) often provide an adequate characterization of the nonequilibrium state. The problem then becomes one of determining the average dynamics of these variables (i.e., macroscopic transport equations) and their fluctuations. At present, almost nothing is known about the transport properties of systems far from equilibrium.
This semiphenomenological approach was first proposed some fifty years ago by Onsager and has essentially been verified for systems near equilibrium.2,3 Recently, it has been extended to apply to
systems far from equilibrium,4,5 but it is difficult to evaluate the success of such an extension since it has been applied to very few problems. The reason it has not been more thoroughly tested is that its application is straightforward only if the nonequilibrium state is characterized by a few simple macroscopic properties.
One such nonequilibrium system is a simple fluid in a state of uniform shear flow. In uniform shear flow, the pressure and temperature are spatially constant while the flow velocity varies linearly with distance in a direction perpendicular to the flow velocity (Figure 1.1). In a rectangular coordinate system (denoted as xyz), the flow field, uo(r), may be written
u (r) = a*r (1.1)
0
where the shear tensor, a, is
a = axy (1.2)
and "a" is the shear rate. The deviation of this system from equilibrium is controlled by a single parameter, the shear rate, giving a simple characterization of the nonequilibrium state.
A sheared fluid is one of the simplest nontrivial nonequilibrium systems imaginable and, indeed, may appear so simple as to be uninteresting. However, even in this simple system drastic effects due to deviations from equilibrium are present, for example in its transport properties, giving the system an intrinsic as well as
Y=L U L U oL
XI
I
Flg. 1.1 Uniform shear flow.
5
pedagogical interest. In the remainder of this introduction, we will first describe three anomalies in the behavior of sheared fluids, and then we will describe the methods and physical ideas which will be used to attempt an understanding of the system. The introduction concludes with an outline of the work to be presented.
Sheared fluids are of value as model nonequilibrium systems not only because of the simplicity of the macroscopic state, but also because a great deal of effort has gone into developing methods of simulating such systems on a computer.613 This is of prime importance because, as will be discussed in detail later, laboratory experiments involving simple atomic fluids are at present limited to the small shear rate regime (i.e., close to equilibrium). As will be discussed in the next chapter, experiments are possible involving more complex rheological fluids, such as polymeric fluids and glasses,14 and colloidal suspensions15 which exhibit many of the properties to be discussed here for simple fluids. We concentrate on simple fluids, however, to simplify the theoretical description. Thus, the only "empirical" data which exist for simple fluids at large shear rates are from simulations. In fact, it is this body of "empirical" results which has spurred most of the recent theoretical interest.4,1625
The first important anomaly observed in these simulations was
associated with the phenomena of shear thinning. As some fluids are sheared, their viscosity drops; they become "thinner" or less viscous. The dependence of viscosity on shear rate had been calculated in 1969 by Yamada and Kawasakil6 as
6
n(a) = n0  lal/2 + O(a). (1.3)
The nonanalytic nature of the dependence on the shear rate was previously unsuspected and indicates the breakdown of nonlinear response based on a formal expansion of the Liouville equation in the shear rate. In addition, the result of the simulations was that n1 was measured to be two orders of magnitude larger than Yamada and Kawasaki predicted. Several other authors16 calculated n1 and obtained different numbers, although all are two orders of magnitude too small. A recent proposal by van Beijeren26 for the resolution of this problem will be discussed below.
In 1982, a second anomaly observed in simulations was reported by Erpenbeck.12 He found that, at very high shear rates, the system passed from the amorphous fluid phase into an ordered phase in which the particles line up one behind another in hexagonally packed "strings" (Figure 1.2). There has recently27 been some debate about the actual nature of the string phase but it is generally agreed that a high shear rate instability, or dynamical phase transition, exists. We shall present below a model exhibiting an instability on atomic scales which matches up qualitatively quite well with that seen in the simulations (note that macroscopic hydrodynamics is believed to be stable for planar shear flow28).
A third anomaly has been reported by Evans7 as occurring in twodimensional fluids. His simulations show the shear viscosity remaining constant as the shear rate increases until at a critical
cr 0
O 0
O
O
00 0 0 a >a*: O O
Fig. 1.2 Large shear rate orderdisorder transition where ac is the
critical shear rate.
8
shear rate the fluid begins to thin in the expected manner. Evans attributes this to a small shear rate instability for which no theoretical understanding exists as yet. However, the transport properties of two dimensional fluids are anomolous even in equilibrium (two is the critical dimension) leading one to believe that Evans' observations may be related to this fact.
Most of the work to be presented here will consist of an attempt to construct and analyze models of fluctuations around the nonequilibrium state in dense sheared fluids. Once the fluctuations are statistically characterized, one is able to calculate time dependent correlation functions which in turn can be related to observable physical properties of the system. The phenomenological ideas behind this approach are very similar to those in the EinsteinSmolochowski model of Brownian motion.29 One considers the dynamics of the hydrodynamic variables (analogous to the position of the Brownian particle) around the steady state. The hydrodynaniic variables are appropriate to describe properties on a scale which is large compared to the microscopic correlation length. Because these variables are the local conserved densities, it is assumed that they decay according to the macroscopic conservation laws (e.g., their time evolution is governed by the NavierStokes equations) with the exception that the equations include stochastic sources which represent the thermal noise responsible for the creation of the fluctuations.4'30 Finally, the stochastic sources are generally taken to be Gaussian white noise with correlations fixed by a fluctuationdissipation theorem.
Having these equations, one in principle solves them for the
hydrodynamic variables as functionals of the stochastic sources and thus may calculate time dependent correlation functions using the known statistical properties of the sources. As alluded to above, these correlation functions can in turn be used to calculate various physical quantities such as the effect of the steady state on transport coefficients and thermodynamic susceptibilities.
While providing a closed and intuitively appealing model of fluctuations, this description is entirely phenomenological and as such represents an uncontrolled approximation. Several authors5,3133 have addressed the question of the formal derivation and validity of fluctuating hydrodynamics and we shall outline their ideas in a later chapter. The result of these studies is that the foundations and validity of this approach are much better understood.
Thus, one approach to the study of nonequilibrium systems is to replace the microscopic description of the system, as characterized by the distribution function and Liouville's equation, with a hydrodynamic description. A hydrodynamic description of a system is obtained as a result of two approximations. The first is that attention is restricted to the local densities of the conserved variables (e.g., mass, momentum and energy in a simple fluid). On length scales which are large compared to the microscopic correlation length (i.e., the mean free path, ï¿½mfp), these variables are expected to decay slowly since their integrals over the whole system do not change at all with time. Thus, on these length scales, these variables are assumed to approximately decouple from the remaining,
10
quickly fluctuating, microscopic degrees of freedom which give rise to dissipation and thermal fluctuation of the slow variables. Concurrent with this, if k is the wavevector (inverse wavelength or gradient in real space) of the conserved densities, it is required that (kimfp) be small. The equations governing the dynamics of the conserved densities may then be expanded in this small parameter with truncation at second order yielding, for a fluid, the usual NavierStokes equations.
In a sufficiently dense fluid, the mean free path is actually smaller than the typical size, a, of the atoms making up the fluid. In this case, the hydrodynamic description will include terms of all order in (ko) while still expanding in (kimfp). The result is a hydrodynamic description of the (dense) fluid which contains information about the structure of the fluid and is valid on atomic (i.e., ka  1) length scales although still restricted to lengths which are large relative to the mean free path (kimfp << 1). One might think, and indeed it was long thought, that atomic scale fluctuations require knowledge of all small scale fluctuations; after all, the reason for considering only conserved densities at long wavelengths is that they decay on macroscopic length and time scales while all other fluctuations decay on much shorter, microscopic, scales. However, while studying models of small length scale fluctuations, de Schepper and Cohen34 discovered that even in this regime the conserved densities decayed much slower than other variables and gave the dominant contribution, on these scales, to the densitydensity correlation function. This indicated that restricting
attention to the hydrodynamic subspace was a valid means of studying fluctuations even on very small length scales. Thus, hydrodynamic models can be used to study even atomic length scale fluctuations provided the mean free path is small enough. In particular, van Beijeren's explanation26 of the discrepancy between theory and simulation results for the size of the shear rate dependent renormalization of the shear viscosity is based on this generalization of fluctuating hydrodynamics to include atomic scale fluctuations in dense fluids.
The work to be presented below consists of modeling hydrodynamic fluctuations as a means of understanding the nonequilibrium state. In Chapter 2, we will discuss experimental and simulational means of studying steady state shear flow and the connection between the two. The third chapter discusses the formal ideas of statistical mechanics used to construct models of nonequilibrium fluctuations and, in particular, to construct fluctuating hydrodynamics at long wavelengths as well as generalized fluctuating hydrodynamics. This is followed by Chapters 4 and 5 analyzing the long wavelength model, the linearized and nonlinear versions respectively, with the aim of deriving various physically observable consequences of the nonequilibrium state. The sixth and seventh chapters contain, respectively, the development and analysis of the small wavelength model and include a discussion of the hydrodynamics and large shear instability occurring on these length scales. The dissertation concludes with a discussion of the results obtained, the completeness of the present description and future directions.
CHAPTER 2
PHENOMENOLOGY AND SIMULATION
In this chapter, we discuss the basic phenomenology, such as length and time scales and strength of shear, which affects the observability of deviations from equilibrium behavior. We then discuss the possibility of performing laboratory experiments which measure nonequilibrium properties. Finally, we discuss the methods of simulating uniform shear flow on a computer.
Phenomenology
In macroscopic hydrodynamics, one typically has three types of
parameters: thermodynamic derivatives such as the sound velocity, co, transport coefficients such as the kinematic viscosity, v0, and the wavevector, k, the magnitude of which is inversely proportional to the length scale of the phenomena one is interested in. The time scale of hydrodynamic dissipation 0s (vk 2) (in equilibrium, momentum fluctuvok t
ations dissipate as e ). In hydrodynamics, one usually uses perturbation theory treating the wavevector as a small parameter. The dimensionless parameter which actually occurs is (v k/c ) thus requiring that v0k << c0 or v0k2 << c k. Physically, this means that the rate of hydrodynamic dissipation is small compared to the rate of sound propagation. This is consistent with the NavierStoke's
12
13
equations governing hydrodynamics since they are themselves the result of an expansion to second order in the parameter (kï¿½mfp), where 2mfp is the mean free path, and in fact c0 /0 z mfp. Typical values of these constants for water at STP35 are ~ 1.8 * 102 cm 0 a
c0  1.5 * 105 cm/s so zmfp = 10 cm. Thus, the applicability of
6 1 1
hydrodynamics requires that k << 8.3 * 10 cm , or ka << 10
8
where a  10 cm is a typical atomic radius.
When a fluid is sheared, a new time scale, the inverse of the
shear rate (denoted here as "a"), is introduced. When the shear rate is as large as the rate of dissipation, we expect that largd deviations from equilibrium behavior will occur. For this reason, we refer to a  v0k2 as the large shear regime. When studying the long wavelength hydrodynamic model, we will require that a < v0k2 thus including the large shear regime while allowing us to use perturbation theory (since, given this bound, a << c0k).
Experiments
Various methods exist for measuring the static and transport
properties of fluids. For example, neutron scattering can be used to measure the static and dynamic structure factors (related to the pair correlation function and transport properties, respectively) on atomic length'scales.36 Light scattering samples the dynamic structure factor on hydrodynamic length and time scales.37 These methods and others are available to measure the corresponding properties for nonequilibrium fluids. However, the present ability to test theories of sheared systems is limited by the shear rates which can be obtained
14
in the laboratory. At present, the largest shear rates reported are less than 104 HZ.38 This shear rate is in the large shear regime
1/2 2 1
when k  (a/vo) ~ 7.45 10 cm . This means that until we look at length scales as large as about 102 cm the effects of shear will be small compared to the usual equilibrium properties. Thus, light scattering, which measures properties on a length scale of the incident beam's wavelength  105 cm, and neutron scattering, which is not practical for wavelengths above about 108 cm, require shear rates about two orders of magnitude larger than those attainable in the laboratory.
There are systems, other than simple fluids, in which similar nonequilibrium effects occur and which have the practical advantage that they have long relaxation times. Because the large shear regime occurs for a  inverse relaxation time ( v0k2 in a simple fluid) it thus becomes accessible in these systems. One example of such a system is a colloidal suspension in which the relaxation time t  d2/D for a typical interparticle spacing, d, and diffusion constant of the spheres, D. In this system, the density of spheres determines the relaxation time and the nature of the probe is not restricted as in a simple fluid. Experimental studies of sheared colloidal suspensions have been performedl5 and numerous nonequilibrium phenomena are observed including a string phase similar to that seen in simulations of simple fluids.
Another class of systems exhibiting easily accessible
nonequilibrium phenomena is macromolecular (e.g., polymeric) fluids
15
and glasses.14 The polymeric systems are not simple fluids because they are composed of anisotropic molecules. However, Simons et al. have shown that strong similarities exist between the phenomena seen in simple and complex fluids. These fluids demonstrate a wide variety of nonNewtonian behaviors including shear thinning. The accessibility of these phenomena again arise due to the fact that the constituents, the macromolecules, have long relaxation times lowering the lower bound of the large shear regime.
In this work, we will restrict our attention to simple fluids.
As the preceding discussion indicated, we will not have direct experimental tests of our results. The reason we do not consider colloidal suspensions, polymeric fluids or some other such system is that the results from computer simulations of simple fluids have shown that even our understanding of the simplest system is incomplete. However, because effects,.such as shear thinning and the string phase, seen in the simulations are observed in physically related laboratory systems, we have some confidence in the simulations. Thus, we propose to attempt to understand the simplest sheared system and hope, in the process, to gain some understanding by analogy of the more complicated systems in which the nonequilibrium state has important physical consequences.
Simulations
Most of the data for simple fluids have been generated by
nonequilibrium molecular dynamics613 rather than in the laboratory.
16
Such simulations require both the establishment of the linear shear profile as well as the maintenance of constant temperature.
In equilibrium molecular dynamics, one considers a box with
periodic boundary conditions containing N particles. At each time step in the simulation, one solves the equations of motion (e.g., Newton's equations) for all of the particles and moves them accordingly. The periodic boundary conditions have the effect of decreasing the size (number) dependence of the results allowing one to achieve the bulk, or thermodynamic, limit with on the order of 103 particles. In this way, one may study static properties, such as the pair correlation function, and dynamic properties like the decay of correlations. Clearly, simulations have the advantage over experiments in that one can obtain much more data about the system since one has the microscopic state at one's disposal., They are, however, limited by the number of particles, N, and the length of time which can be simulated. In particular, the size restriction limits the smallest wavevectors which can be studied.
The linear shear profile is generated by modifying the simple
periodic boundary conditions of equilibrium simulations. The system remains periodic in the x and z directions but when a particle with momentum p passes through the y = L plane at time t and with coordinates (x,L,z), it is reentered at y = 0 with coordinates (xaLt, y, z) and momentum 2  maLx. These boundary conditions are known as LeeEdwards' boundary conditions. They become more transparent when expressed in the local rest frame coordinates, q and p, defined as
17
. = 2  a t
. = 2  ma . (2.1)
In these coordinates, the local flow velocity is zero and it is seen that the LeeEdwards' conditions correspond to simple periodic boundary conditions in the rest frame.
The temperature, T, is a measure of energy fluctuations around the steady state and is defined by
3 1 N
NkBT = < (i  ma. * q)2> i=1
N
<1 (Pi)2>, (2.2)
where kB is Boltzmann's constant. Substitution of the conditions of spatially uniform shear flow into the NavierStokes equations shows that a sheared fluid will undergo viscous heating. Specifically, one finds that
_E 2
a= na , (2.3)
where e is the internal energy density, leading to
E(t) = E(O) + na2t. (2.4)
18
For water sheared at 104 Hz at STP this leads to a heating rate of 0.024 K/s. For temporally uniform shear flow, particularly for the very high shears applied in computer simulations, an external, nonconservative force is necessary to hold the temperature constant. This heat extracting interaction is called a thermostat. This is equivalent to adding a drag force to the equations of motion:
d
d
dt = i x(r) i
2
(r) = J (F.  a ) / 2. (2.5)
J J
where the form of A(M) has been chosen to hold the temperature as defined in eq. (2.2) constant.
The equation of motion form for the thermostat is important
because, with the LeeEdwards' boundary conditions for imposing the flow field, the simulations can be mapped onto a well defined mathematical problem: the solution of Liouville's equation, modified to account for the thermostat, with the imposed boundary conditions. The modified Liouville equation is obtained from the general form of Liouville's theorem,39
a 1 a a a a + + F * ]p(r,t) ai [A(r)ip(r;t)]= 0
Sq, a(2.6)
(2.6)
19
where the repeated index, i, is to be summed and runs from 1 through N and the dot indicates an ordinary vector scalar product. So, (2.6) may be rewritten as
 p + Ep = 0, (2.7) with
A A
ï¿½ = + fint + zext (2.8) where the free streaming and interaction terms are, respectively,
0 i 3q,
ï¿½ = F i (2.9) int i Dp
i
while the external term, due to the nonconservative thermostat, acts on the distribution function as
ext p(r;t) = [(r)i1p(r,t)]. (2.10)
Because this operator includes the nonconservation of volume in phase space due to the nonconservative force (Liouville's theorem states that d (p(r)dr) = 0 and with nonconservative forces  dr * 0), it has a different form when acting on any other phase function.
20
Specifically, if A(r) is an arbitrary phase function, we have that
ï¿½ A(r) =  ) * A(r). (2.11)
exti ï¿½ i
The bar indicates the difference between these operators and we note that ï¿½0 = ï¿½0. It is easy to see that ï¿½ and ï¿½ are adjoints with respect to phase space integration.
Finally, we should mention that when working with hard spheres, the potential is singular and the interparticle force, Fi, appearing in these equations ill defined. This difficulty is eliminated by replacing the interaction term with the socalled Toperators which generate the singular hard sphere dynamics. Henceforth, for hard spheres, we will write
^+
ï¿½int  Tï¿½(,pi'j;qij) (2.12)
where ï¿½ refers to forward and backward time equations. The form of these operators can be found in the literature40 and, for example,
T(x1,x2)A(1 ,x2) =  a2fd; (( o * ~12 12 12
* [b12  1] A(x,X 2)} (2.13)
where x, = a Iq, etc.; a is the hard sphere diameter, 12 =  2
21
and
b AA AA
bl2A( l' 1; 2'2) A(ql' 1  o ~ ;l2 '2 + " 12). (2.14)
We note here that the Toperators are not self adjoint, as is the interaction term for continuous potentials, eq. (2.9), so their adjoints are denoted with bars as well.
CHAPTER 3
FLUCTUATION MODELS
In this chapter, we discuss the statistical mechanical basis of the models to be analyzed. Specifically, the models we will study will be seen to be two different limits of a general FokkerPlanck equation. The detailed form of the models will also be given.
Generalized Langevin Equation
Our basic approach, as discussed in the introduction, to studying nonequilibrium fluctuations around shear flow is to choose a restricted set of variables whose dynamics we wish to study. At a fundamental level, any classical system of N structureless particles, or atoms, is described by the 6N coordinates, collectively called r, which are the positions, qi, and momenta P, of the atoms (where i ranges from 1 to N). The distribution function of the system, p(r;t), which describes the density of states in phase space, satisfies the Liouville equation,
t p(r;t) + ï¿½ p(r;t) = 0; (3.1)
the explicit form of which was discussed in a Chapter 2. The set of restricted variables we choose to study is the set of conserved
22
23
A A
densities. These are mass density, p(r), energy density, e(r), and momentum density, j (r). A caret over a function indicates that it is a phase function. Their definitions are
S N
p(r) = Z m A(r_)
n=1
^ N N
e(r) = E  _ A(rqn) + F Z V( I,j I) A(r _)
n=1 nnn'
A N
ja(r)  Z pn A([n) (3.2)
n=1
where A(r) is some localized function of r (for instance, a Dirac delta function). These variables will be denoted collectively as y(r). Our reasons for choosing this particular set were alluded to in the introduction. At length scales large with respect to the mean free path, we expect that the conserved densities are the only variables with long relaxation times (i.e., long with respect to the collison frequency). Therefore, we expect that the statistical properties of the system at these scales can be adequately described in terms of these variables with the other degrees of freedom being modelled as rapidly varying sources. In dense fluids, the mean free path can be smaller than the atomic dimensions and it is found, as one might guess, that the hydrodynamic variables, or conserved densities, still have long relaxation times and so still, on the average, obey an approximately closed dynamics. Because the other degrees of freedom are more rapidly varying, they may again be modelled as stochastic sources.
24
In the analysis to follow, it will prove convenient to work with the Fourier transform of the variables in eq. (3.2). This is defined to be
^ =ikfr
1
y(k) = VI dr e y(r). (3.3)
We also note that the vector y(k) is explicitly,
y(k) = (p(k), e(k), (k)). (3.4)
For notational convenience, we will write
ya + y (k) (3.5) where the Latin index includes both the vector index a and the wavevector k. Repeated Latin indices are to be summed. The average of ya, or y (k), over the nonequilibrium ensemble, p(r;t), will be denoted as where
f J dr p(r;t)ya (3.6) and are the macroscopic variables one would measure in the laboratory.
In general, one also defines time dependent phase functions, y(k;t), as
25
ï¿½t
ya (k;t) = e y (k) Ya(t) (3.7) where
= fdr Ya e p(r;0) = fdr [et ya] p(r;O) = fdr Ya (t) p(r;0) = . (3.8) So, the macroscopic variables are the average over an initial ensemble, p(r;O), of a (t). The fluctuations of a (t), denoted as 6ya(t), are defined to be 6Ya(t) = Ya(t)  . (3.9) Because we are restricting attention to the set of variables y(t), it is useful to introduce the distribution function of these variables, p(y;t), defined to be
26
p(y;t) = fdr' 6(yy') p(F';t), (3.10) where y is some possible value of the phase function y. So, for any function of the y's, say A[y], it is clear that
= fdr A[y] p(r;t)
= fdy A[y] p(y;t). (3.11) Finally, for a steady state such as uniform shear flow, we have that pss(r;t) = pss(r;0)  Pss(r). So, (y a;t> = ss
= Yo0a (3.12)
Grabert et al.41 show, using the projection operator technique introduced by Zwanzig,42 that the distribution p(y;t) satisfies the following exact generalized FokkerPlanck equation, a p(y;t) + a v [Y] P(Y;t)  fvads fdy' Daby,y';s] ps(Y') aya
aay
p(y';ts) } = 0 (3.13) SYb' P S(y')
27
where pss(y) is the steady state distribution, Va[y] is called the drift vector and D aby,y';s] is called the diffusion matrix. Repeated Latin indices are again to be summed in (3.13). The functionals v and D are expressed in terms of steady state averages and their detailed form is given in Grabert et al.41 Here, we only note that
 {Va [y] p (y)} = 0 (3.14) ayaa
so p s(y) explicitly satisfies (3.13).
Associated with (3.13) is an exact generalized Langevin equation for the variables Ya(t),
t Ya(t) a  aby(t)] + ds DabY(ts);s] Fb[Y(ts)]
a A
 A Dab[y(ts);s]}
ayb(ts)
= R (t) (3.15)
where the thermodynamic forces, Fb[y(ts)] are defined in terms of the steady state distribution as
SIn p s(Y)
Fb[y(ts)] =  = y(ts (3.16)
which also appears in (3.13). The functional va is the same as that appearing in (3.13) while
28
Daby;s] = fdy' D [y,y';s] l (3.17)
ab ab y=y.
The sources Ra(t) have the following constrained averages,
0 = fdr' 6(y y ') R '(t) p (r ')
D b[Y;t] = dP' 6(y y ') Ra'(t) Rb'() p (r') (3.18)
where the hat on Ra(t) indicates it to be a phase function and where a primed function has argument r'. The first of equations (3.18) indicates that Ra,(t) is orthogonal, under the steady state average,
A
to any function of y and indicates that the constrained average of R (t) at fixed y is zero. The second of equations (3.18) has the form of a fluctuationdissipation theorem relating the correlation of the "fluctuating force" to the diffusion matrix.
Now, eq. (3.15) together with properties (3.18) may be modelled as a stochastic differential equation if the sources R (t) are idealized as stochastic processes, denoted as Ra(t), which have the form,
Ra(t)  Mmnly] Sn(t) (3.19) where
=
= (27) 6(k+k') 6(tt') A
A aa 6(tt') (3.20)
29
where A,' is a constant matrix and the amplitudes Mmn[y] are fixed by the fluctuationdissipation theorem, eq. (3.18). Finally, the variables Ya(t) are replaced by stochastic variables denoted as a (t) = y (k,t) and the average over the distribution function replaced with that over the stochastic sources,
ss a R(t) (3.21)
since the variables Ya(t) may be considered to be functionals of the stochastic sources.
Long Wavelength Fluctuations
If we restrict our attention to long wavelength fluctuations, then the preceding Langevin model can be explicitly evaluated. Specifically, we only consider wavevectors such that
kimfp << 1. (3.22)
In this case, the various functionals appearing in eq. (3.15) may be evaluated in terms of a power series in ktmfp. If we truncate the expansions at second order, then the model may be given explicitly,5 in terms of the usual microscopic conservation laws,
atp + V 0
30
at
Sj + t = 0 (3.23)
where sa and taB are the heat flux vector and stress tensor, respectively. These in turn each have a convective part, a dissipative part and a stochastic part, where the latter represents the stochastic source,
8 R
s = eu + t uB + A T+s a = a a o 3r a
t = u ( u )  V + t (3.24) ta =uJ + P  aB Dar o aB a
where we have introduced the fluctuating velocity, ua = p j , the fluctuating pressure and temperature, p and T, and transport coefficients, n , Ko and A . These fluctuating thermodynamic functions and transport coefficients are the same functions of y as one has in equilibrium; i.e.,
p[y] = pey]
= pe[p,'] (3.25)
where pe is the equilibrium pressure functions and the last line
31
indicates that pe depends only on the (local) density and internal energy density, E, where E is given by
1 1 .2
E = e  .p (3.26)
The dissipative parts of the heat flux vector and stress tensor are seen to be given by a fluctuating version of Fourier's law and Newton's viscosity law, respectively, where
2
S6 6 + 6 6  6 6 (3.27) aBPv ap av av Sp 3 aB Pv
Finally, as previously mentioned, the stochastic sources are modelled as Gaussian white noise with covariances determined by the fluctuationdissipation theorem, eq. (3.18),
R 2 1/2
S = (k A T ) s (r,t)
a Bo a
R 1/2 1/2
tB = (kBTno0) s1s(rt) + (kBTio) /s2s (r,t) (3.28)
= = = 0
= 2 6 (r r 2) (tlt2)
a 1't)s (r2't2)> = 2 6a 12 1 2
= 2 A 6(r r2) (t t2)
= 2 6 6 6(r r2) 6(tlt2) (3.29)
 1 2 p 26B
32
with all other cumulants equal to zero. Equations (3.23)  (3.29) completely specify the fluctuating hydrodynamics model which may be summarized as a Langevin equation for the conserved densities, Ya(t), with a deterministic part given by the NavierStokes equations, local equilibrium transport coefficients and susceptibilities supplemented by the stochastic sources, eqs. (3.28) and (3.29), which model the excluded degrees of freedom. Finally, we note that this model is given in the language of eq. (3.15) in Appendix A where it should be noted that
D aby;t] = 0 (3.30) ayb
as is shown by Zubarev and Morozov.5
Small Length Scale Fluctuations
In a dense fluid, the mean free path may actually be smaller than the typical atomic length scale, to be denoted here as a (for a gas of hard spheres, a would be the hard sphere diameter). In this case, one may study atomic length scale fluctuations with the present methods. However, the evaluation of eq. (3.15) requires care for while we still expand in (k mfp) we must keep all orders in (ka). Quantities dependent on (ka) carry information about the small scale structure of the fluid. We can evaluate (3.15) while retaining the desired information regarding the small scale structure by taking its short time limit to obtain
33
at ya(t)  Va[Y(t)] = Ra(t). (3.31)
The drift vector is defined to be
Va[Y(t)] =  {<6(yy) ï¿½ Ya>ss/<6(yy)> ss}i . (3.32)
In general, this is a complicated nonlinear functional of y(t). However, as is explained in Chapter 5, we will only be interested in this functional expanded to second order in the deviations from the stationary state. Thus, we expand as
v [y(t)] = v[y (k)] + 6Yb(t) * [ A v [y(t)]] + etc.
6yb(t) 0
(3.33)
The first order term is
6 1
S v a[y(t)] = (<6y6y> ) b<6y y > ss (3.34)
y b(t)
This model is explicitly evaluated in Chapter 6 and analyzed in Chapter 7. While it may seem strange to study the short time limit of the generalized Langevin equation, it should be noted that similar short time equations have been found to provide qualitatively accurate descriptions of dense fluids even out of the short time regime (an excellent example is the Enskog equation43).
CHAPTER 4
LINEARIZED HYDRODYNAMIC MODEL
In this chapter, we linearize the NavierStokesLangevin
equations described in Chapter 3 and use them to calculate various statistical properties of the system. We will do this in such a manner that the results apply for shear rates up to and including the large shear regime as previously defined. Of particular interest will be the hydrodynamic propagators and the two time correlation functions. We will discuss the applications of our results to light scattering experiments and will make contact with the small shear results of other authors. The chapter concludes with a description of various formal statistical properties of the system highlighting the difference of this model from local equilibrium.
Linearized Model
The complete hydrodynamic model was given in Chapter 3, eqs. (3.23)  (3.29), and will not be repeated here. The variables we choose to work with are the mass density, total energy density and momentum. So,
y + (p,e,j) (4.1)
34
35
and we denote the steady state averages of these quantities with subscript naught,
y (P0,'e'0,) (4.2)
and the deviations from the steady state, the fluctuations, are denoted by Za,
Z (r,t)  y a(r,t)  y (r,o) (4.3)
where we also use internal energy density and velocity in place of total energy density and momentum density. Finally, we scale and Fourier transform Z (r,t), to obtain the variable Z(k,t) in terms of which we wish to work,
Za(k,t) = fd3r k*r Z (r,t) (4.4)
Z a(k,t) + c 6p(k,t), 6E(k,t), e (k)u(kt)] (4.5)
2 2 O 2 o 0
c1 = p( ) ; c = (  ) (4.6)
1 0 8p0 0 3C0 0
where p0 = PLE POO] is the steady state, macroscopic pressure; h0 = O + P is the enthalpy density and (e (k)} are a set of ^(1) ^
pairwise orthonormal vectors with e (k) = k/k = k and E is the internal energy density and u = pj is the velocity. The scaling, internal energy density and u = p j is the velocity. The scaling,
36
(4.6), will simplify the form of the equations below and we emphasize that cl and c2 are equilibrium thermodynamic susceptibilities.
In the course of Fourier transforming and linearizing the
equations, we find that we must transform the convective terms which '1
include the operator p0O * V which on an arbitrary function f(r) looks like
1
P00 * V f(r) = a irj  f(r). (4.7)
The transform of which is
3 ik*r 3 ikr r
fdr ik a..r f(r) = a (i )fd3re f(r)
 ajr la k. r
= a k f(k). (4.8) ij 8k i
Thus, it turns out that the form of the linearized equations is
a t
 Z (k,t)  a k Z (k,t) + H (k;a) Z (k,t) = R (k,t) (4.9) Bt a 3 j a  a  a
where Ra represents the stochastic sources which, at linear order, do not depend on the fluctuations, Za. The details of the linearization including the explicit forms of H a and Ra are given in Appendix A. Equation (4.9) is the primary result of this section. In the following sections, we will analyze the hydrodynamic (i.e., normal) modes described by this equation. For zero shear, H is the usual
37
matrix of linearized hydrodynamics the eigenvalues of which describe sound propagation, heat diffusion and momentum diffusion.
An important feature of eq. (4.9) is the coupling of variables of different wavevectors implicit in the gradient term. This unusual feature has important consequences on the evolution of fluctuations and is known as linear mode coupling.
Green's Function
To solve equation (4.9), we introduce the Green's function GaB(k,k';t) defined to satisfy
" G (k,k';t)  a ki G (k,k';t) + H o(k;a)G B(k,k';t) = 0
a i k i a(kaku;t)
with Ga(k,k';O) = (2w)3 6(k  k') 6a. (4.10)
In equilibrium, i.e., a = 0, the solution of (4.10) requires nothing more than the diagonalization of the hydrodynamic matrix H . Before solving (4.10) for nonzero shear, we recall the equilibrium solution.39
Let the eigenvectors of HaB(k,a=0) be denoted as
((i)(k)}ji=1,5 and the corresponding eigenvalues be {A(i)(k) i=1,5 so that
H a (k,O) o (k) = i) (k) i)(k). (4.11) a  a
38
Furthermore, as Ha is in general nonhermitian, we introduce the biorthogonal set of vectors, {In(k)} i=1,5' defined to satisfy
(i) (j)
n (k) . ( (k) = 6... (4.12)
In terms of these functions, the Green's function can be immediately determined to be
G (k,k;t) (k) )(k')(2) (kk')exp  [Ai(k)(tt0)] i=1
(4.13)
and tO is a constant which can be fixed by an additional boundary condition. The solution to eq. (4.9) is then given by
dk'
Z(k,t) = f 3 G (k,k';t) Z (k',t )
(2r)
dk'
+ dt d f ) G (k,k';tT) R (k',) (4.14) t 3(2n)3 as 8
where the wavevector integrations extend to infinity.
Because the hydrodynamic matrix H.0 is the result of a truncation at second order of an expansion in (kimfp), the eigenvalues can be determined by perturbation theory in this same small parameter. In equilibrium, one obtains two sound (propagating) modes, a heat mode and two degenerate shear modes,
39
1 2
A (k) = ic k + 2 rOk
1 2
A2(k) = ic k + 2 r k
A (k) = DOT k
A4(k) = VOk2
A 5() = 0 k2 (4.15)
where explicit forms for the speed of sound, co, the sound damping constant, r0, and the thermal diffusivity, DOT, will be given later. For the present, we note only that the modes all decay like exp(constant)k2t and that these results can be used to calculate the intensity of scattered light from an equilibrium fluid (as described below). Upon doing so, one obtains the classic LandauPlaczek formula44 which is perhaps the most successful application of fluctuating hydrodynamics. In fact, comparison of this result and experiment provides one of the best known means of measuring the equilibrium transport coefficients.
Returning to the nonequilibrium problem, eq. (4.10), we introduce a shear rate dependent generalization of eqs. (4.11)  (4.12) as follows:
40
a k  ( (k;a) + H (k;a)( (k;a) = M (k;a)( (k;a) ijk ak a  a a
(i) (j)
n (k;a) * (k;a) 6 . (4.16) 1j
The presence of the gradient operator in (4.16) greatly complicates the problem. However, as discussed in Chapter 2, we restrict the size of the shear rate by requiring that the rate of convection be no larger than the rate of dissipation yielding the condition,
a < vk2. (4.17)
This choice is what allowed us to neglect heating in these equations,
2 14
since the rate of heating goes like a  k (see eq. 2.3), and makes possible a perturbative solution of eq. (4.16). The details of this calculation are given in Appendix B. The eigenvalues so obtained are
1 2 2 A = ic k +  (r k + ak k /k )
1 0 2 0 x y
1 2 2
A = ic k +  ( k + ak k /k 2)
2 0 2 0 x xy
2
A3 = DOT k
4 = vgk2  akxk /k2
A5 =~vOk2. (4.18)
41
We see that the degeneracy of the shear modes has been lifted while the sound damping constant has been shifted and the heat mode remains unaffected at this order.
If we now write the Green's function as
G (k,k';t) = 1 n (k;a)( (k';a)G )(k,k';t), (4.19) i= a
then we find that the function G(1) satisfies the following equation:
 a jk k + A (k,a)] G (k,k';t) = 0
Bt ij i ak.
i
G (k,k';0) = (2w) 6(kk'). (4.20) This scalar equation is easily solved and if we introduce a time dependent wavevector as
k.(t) = k.  kjaji t = k.  6iy k xat, (4.21) then the solution of (4.20) is
G i)(k,k';t) (27)3 6(kk'(t))E(i)(k,t)
where E (i)(k,t) = exp  t d A(k(T)). (4.22) The explicit form of the propagators in (4.22) is
42
E(1)(k,t) = [k/k(t)]1/2 exp{ic ka(t)  1 2 )
S(1 ) 1/2 .k/kt 2 0 t
E (2)(k,t) = E (k,t)
E(3)(k,t) = exp{DOT k 2(t)
E (k,t)  [k(t)/k] expk(t)
E(5)(k,t) = expfvk 2B(t)} (4.23) where a(t) = (2ak xk)1(k y(t)k(t)  k yk) k(t) + k (t) sgn(ak ) k + ky sgn(akx) y x
8(t) = (akk2 )1 {k2(k (t)  k) + 1 (k3(t)  k)
2 2 2
k = k + k . (4.24)
x z
These functions behave for long times as a(t) + (ak /2k)t2
1 23
8(t) + 1 (ak /k) t (4.25)
3 x
43
Because lim B(t) > 0, the modes are all stable as expected (see Drazin and Reid and Case28 for further discussion of the stability of
macroscopic shear flow).
In order to illustrate the physical significance of some of the shear rate dependent effects described by these equations, consider the evolution of a purely hydrodynamic sound excitation along the direction of flow. Specifically, we take Z (k,) () 6(kk with k = k0x. The evolution of hydrodynamic excitations is given by eq. (4.14) when conditionally averaged over the stochastic force with fixed initial condition Z (k,0). Denoting this average by a double bracket we have that
dk'
<> = f 3 Ga (k,k';t) Z (k',O)
(21)
= 6(k  k (t)) (1)(k) E ()(k,t). (4.26)
Inverting the Fourier transforms yields,
(1) 2 1/4 Ok8(t,a)/2 <> = E (k0(t))(1 + a2t2)/4 er0B(ta)/2
x cos{k [x(t)  0 a(t,a)]} (4.27)
where a(t,a) and B(t,a) were defined in (4.25) and x(t) = x + ayt. The time dependent coordinate x(t) is due to the expected Doppler shift in the sound wave which travels in a moving medium. This kinematic effect may be suppressed by restricting attention to the
44
y = 0 plane and without loss of generality, we may chose x = 0 as well. Thus, the remaining shear effects are due to gradients of the velocity field and not the magnitude of the velocity. Figure (4.1) shows <>/<> for the choice of parameters (ck0/Okr2) = 10 , which is typical, and the two cases (a/r kO) = 0 and 1. Note that in the former case,
rOk2t/2
<>/<> = e cos(c0 k 0t) (4.28)
which is the expression for ordinary sound propagation. Two effects are easily seen. First, the effective sound velocity increases with time relative to that for zero shear rate on this time scale. Specifically,
Ck a(t) 1 21/2
c k {( + T2)/2+ ln[(1 + T ) + ]} (4.29)
C00kot 2
where T R r 0k t. Second is that the attenuation is enhanced, relative to that for zero shear, as time increases,
r0k08(t,a) 1 2 2
= 1 + a . (4.30)
0 0
These relative magnitudes of convection and dissipation would be very difficult to reach with simple fluids but are possible for laboratory colloidal suspensions.
Fig. 4.1 Propagatiqn and damping of sound wave in shear flow for a/rFkO = 1 () and
for a/rk0 = 0 ().
46

.
e un
qllm tl ' ï¿½ n l
milIII a ilirmt lb l! ii 411o Iiib
r'al ~l
47
Correlation Functions
The correlation functions of the deviations of the fluctuations from the steady state are useful for a number of reasons. From them one may calculate the dynamic structure factor, which describes the intensity of light scattered from the system and transport properties such as the shear rate dependent renormalized viscosities and they characterize the deviation of the system from local equilibrium. In this section, we calculate the two point correlation functions.
The two point correlation functions are defined to be
C a (k,t'k',t';a) = (4.31)
where the brackets denote the nonequilibrium average which in the present case corresponds to an average over the stochastic forces. Time reversal invariance and stationarity require that
C s(k,t;k',t';a) = Cas(k,tt';k,O;a)
C a(k,t;k',t';a) = P P C (k',t;k,t'a) (4.32)
where Pa = ï¿½1 is the parity of Za under time reversal. Therefore, we may consider t > 0 without loss of generality.
If we define
(k,t) = n W(k;a) R (k,t), (4.33) a a
48
then, using eq. (4.14), we find that Ca (k,t;k',O;a) Ca (k,k';t) = fO dTfodT{ (i) (k)5 )(k')E (k,T)E ()(k',T')
i,j
x (4.34) where we have taken tO +  in eq. (4.14) and used the boundary condition that Z (k,) is finite. It is clear that
= (2) 36(k+k')6(tt')F (i)(k) (4.35) where the form of F(i) (k) is easily determined. Thus,
3 (1) (i) (i) C (k,t;k') = (21) 6(k+k'(t)) (i) (k)n (k')E (k,t)C (k';a).
(4.36)
The functions C S(k;a) are related to the equal time condition functions because, in particular,
CaB(k,t=O;k';a) = (2)36(k+k')C a(k;a). (4.37) The explicit form of these functions is
49
C (k;a) = Ir ( { (i) E(j (k)E() (k,)E() (k,T)
ae 0 a 8
i,j
x F(ij)(k(T))}. (4.38)
The nonzero elements of C a(k;a) are given in Tables 4.1 and 4.2.
Some immediate information may be gained from these
expressions. In particular, eq. (4.36) shows that the two body correlations decay in the same way, i.e., with the same propagator E(i)(k,t), as the macroscopic modes. This represents a generalization of Onsager's regression hypothesis to nonequilibrium states. Furthermore, eq. (4.38) is a fluctuationdissipation theorem for the nonequilibrium system. These statements are more obvious when cast in the form of differential equations:
[ a k, ]C (k,M';t) + H (k;a) C (k,k',t) = 0 (4.39)
at ii a1j k  o aS
[6 a.k.  + H (k;a)] C (k;a) + [H (k;a)] C (k,a) = R (k) aa ij i 3k. a aS  aa ao a
3
(4.40)
where RaB(k) is the amplitude of the fluctuating forces,
s 6(k+k') 6(tt') R (k). (4.41)
Equation (4.39) shows more clearly the linear regression law while (4.40) has the usual form of a fluctuationdissipation relation.
50
Table 4.1 Equal time correlations.
2
C (k;a) = kB x [1 + A (k,a)] pp B 0 0 T 1C p(k;a) = kBTOhOOx T[1 h T T + A1 (k;a)]
1 T2 2 C (k;a) = k T x Ex T oC + [h  T ] + h A (k;a)] CE B 0 T T 0 v 0 0 xT 01
1
C44(k;a) = kBToP0 [1  A2(k;a)] C45(k;a) = k BTOPO A3(k;a)
(
C55(k;a) = k T0P0 [1  A (k;a)] 55 B 0 0 4
51
Table 4.2 Nonequilibrium contributions to equal time correlations.
kk k (t) r0 k (t) A(k;a) = aY fdt e k3 (t)
k k (t) 2 k(t) A2(k;a) = 2a odt X  e
k
k k (t) k 2v0k2 (t) A (k;a) = aodt[2xky F(kt) +
A4(k;a) = 2aodt F(k,t)[ x y t F(k,t) + le
k(t)
F(k,t) = A(k(t)) A(k)
kk k
A(k) = z tan 1 Y)
k k k
x
k(t) = k  k * a t
2 2 2 k2 . k + k
x z
52
Physical Properties
We now discuss two physical properties which may be determined
from the correlation functions. They are the long range nature of the correlations and the light scattering spectrum.
The densitydensity equal time correlation function in real space is easily seen to have the form
dk
 ik*r
C (r;a) = f ) e   C(k;a)
pp (2)
= 3kB TO PXT[6(r/t) + /r F(r/i)] (4.42)
where the isothermal compressibility, XT = p ( )T' and where I = (r /a)1/2 is a typical length scale associated with the shear rate. The limit F(r/t)=O yields the usual equilibrium result where the 6function is an idealization of the atomic scale correlations which actually exist in the fluid. The part porportional to r/k has been obtained before via an expansion to first order in the shear rate.20 Because of the form of F, this actually requires that r/k << 1 and yields
F(0) = a rirj /8rYar2 (4.43)
where Y = cp/cv is the ratio of specific heats. This does not, as has been claimed,20 however, allow us to conclude that the correlations are longranged because it is the product of a short distance expansion. To determine the true nature of the correlations, at large
53
distances, we consider a somewhat simpler problem than the larger expansion of F(r,t). Specifically, we consider
fdyjdz C (r) T kBTOPOXT[C(x/) + YH(x/t)]. (4.44)
The part of the correlation function due to the deviation from equilibrium, H(x/t), is given by
H(x/1)  (2w)1 fO dt{(1 + t3/2 /2 (1 + t2)
0 3
x exp[(x/i) (4t(1 + t2)) ] (4.45)
and for large (x/) behaves as
H(x/i) 4 0.828(x/)5/3. (4.46)
At short distances, this function is Gaussian and the complete function is shown in Fig. 4.2. It is thus clear that the correlations are in fact long ranged, decaying as a power law.
The second physical property we wish to study with the
correlation functions is the intensity of light scattered from the fluid. It can be shown39 that the intensity of light at the point r with frequency w due to the initial scattering of a coherent beam of wavevector kO and frequency w0 is given by
0.6
I
H(X/1)
0.4
0.2
0 2.0 4.0 6.0 80 10.0 X/i
Fig. 4.2 Nonequilibrium contribution, H(X/1), to densitydensity correlation function
along direction of flow for a/rOkO = 1 () and asymptotic form ().
0 0o
55
1 2 2
I(r,w) = j 10 w (a/2 C0r)sin 0 S(k,_) (4.47)
where I0 is the spectral intensity of the incoming beam, a is the polarizability of the particles composing the medium, r * EO a EO cos, ED is the amplitude of the coherent beam, k = kj  (W/CO)r and Q = w  O0 and where the structure factor, S(k,Q), is defined as
dk dk
f'*dt d1 2 int S(k,Q) = fdt (kk )8(k+k )C (k t;k ,0;a)e . (4.48)
6 1 2 pp 1
* (21)
Thus, we see that the densitydensity correlation function determines the intensity of scattered light. The functions, j(k), are Fourier transforms of form factors which limit the volume of integration to that irradiated. In Appendix C, it is shown that S(k,0) can be written in terms of nonequilibrium Brillouin and Rayleigh peaks:
S(k;Q) = S R(k;Q) + S B(k;Q) + S B(k;) (4.49)
where
1 a 1 2 (1) (1)
S (k;a) = 2Re{[it + X (k;a)   ajk ] C n (C (k;a)l B 1 2 ij k 1 1 a al
1 a 1 2 (3) (3)
SR(k;Q) = 2Re{[it + X3(k;a)  1 ai k  C1 na C (k;a)}.
(4.50)
56
In equilibrium, the gradients are not present and these
expressions reduce to simple Lorentzians. Expansion of eqs. (4.19) to first order in the shear rate yields previously obtained results; however, the general evaluation of (4.19) is difficult. A quantity which can be calculated directly is the LandauPlaczek ratio, R, defined as
R(k;a) i dw SR(kw)/fdw[SB(k,w) + SB(k,w) S[Cpp (k;a)  Polh0 Ccp(k;a)]/a2C pp(k;a) + p0/h0 C p(k;a)]
(4.51)
where a = p0c1 /ho02. In particular, the integrated Rayleigh intensity is unchanged from equilibrium, fdw S (k,w;a) = (Y  1)/Y (4.52) where Y = ratio of specific heats. An explicit form of (4.51), obtained using the expressions for the equal time correlation functions given in Appendix C, is
R(k,a) = (Y  1)/[1 + YA (k;a)] (4.53)
57
where A = [C (k;a)  C (k;0)]/C (k,0). In Fig. 4.3, we plot
1 pp  pp  pp
AR  [R(k;a) _R(k;O)]/R(k;O) as a function of e for
2
k = k (cosex + siney) and a/r Ok = 1. The small shear rate expansion of these results is in agreement with earlier work by Kirkpatrick et al.20 The result, eq. (4.53), is unexpected in that it states that the integrated intensity of the Rayleigh peak is unaffected by the shear.
Probability Densities
In this section, we recast our previous results in the language of the probability distributions of the fluctuations. In particular, the probability distributions are more useful than the description in terms of fluctuations for the calculation of higher order correlation functions or for averages of more general functions.
The probability and joint probability density are defined by
P(Z,t) = <6(Z(t)  Z)>
P(Z,t;Z',t') = <6(Z(t)  Z) 6(Z'(t')  Z')> (4.54)
where the averages are taken over the stochastic forces, {Ral, and over a specified ensemble of initial values, {Z(0)}. The 6functions are short for
6(Z(t)  Z) H 6(Z (k,t)  Z (k)). (4.55) a,k
0.30
0.10
U /4 I 34 n
Fig. 4.3 Deviation of the LandanPlaczek ratio, eq. (4.52), as a function of angle
for a/nk 2= 1 and k = 0. The value at e = /2 corresponds to k = 0 for which the
0 z x
equilibrium result holds.
CO
59
We begin the evaluation of eq. (4.54) by calculating a related quantitythe conditional probability density defined by
W(Z,t;Z',t') = P(Z,t;Z',t')/P(Z',t'). (4.56) Because the system is being modelled by a Markovian process, the conditional distribution is independent of the initial values, {Za ()} and because of the neglect of heating, the system is stationary so W(Z,t;Z',t') = W(Z,tt';Z',O). So, without loss of generality, we consider
W(Z,t';Z',O) = <6(Z(tIZ')  Z>R (4.57) where, as indicated, the average is only over the stochastic forces and Z(OIZ') . Z'. From eq. (4.14),
(tdZ') = G(t)Z' + dT G(tT) R(T) (4.58) where we have suppressed internal summations. An integral representation of the 6function in (4.57) and use of (4.58) yields W(Z,t;Z',O) = fdA eiA[ZG(t)Z']
 fdA exp{iX[Z  G(t)Z']  12M}
60
= Idet[M(t)]1 i2exp{ [ZG(t)Z'] * M1(t) * [ZG(t)Z']} (4.59)
where the second line follows from the first upon use of the Gaussian statistics of the source term. The matrix M(t) is given by
M(t) = fdT dr' G(tT) * * GT(tT'). (4.60) Equations (4.39) and (4.40) (the regression law and fluctuationdissipation theorem) allow us to reexpress (4.60) in terms of the correlation function matrix, C(t) C ag(k,t;k',0),
M(t)  C(0)  C(t) C1 (0) * CT(t). (4.61)
We can now immediately write down the probability and joint probability densities,
P(Z,t) = fdZ' W(Z,t;Z',0) P(Z',0)
P(Z,t'Z',t') = W(Z,t;Z',t') P(Z',t'), (4.62) where P(Z',O) is the initial probability density. From (4.61), we can determine the stationary state,
61
P (Z) = lim P(Z,t)
ss
t**
= {det[rC(O)]}1/2 exp[ Z C 1(0) * Z (4.63)
which follows from the limits,
lim C(t) = 0 => lim M(t) = C(O) (4.64)
t** t**
and, lim G(t) = 0. Thus, as expected, because the driving process was
t.*
Gaussian and the equations linear, the distribution of fluctuations is also Gaussian. Furthermore, eq. (4.63) shows that the same asymptotic stationary state is reached regardless of the initial condition.
All of these results could have been obtained if we had worked with the FokkerPlanck equation rather than the fluctuation equation. For completeness, we write down the FokkerPlanck equation for P(Z,t) which is easily obtained. For example, by differentiating (4.62) with respect to time,
 P(Zt) J az J (Z't) + J1(Zt)] (4.65)
where the drift term is
J (Z,t) = [a k.  Z (k)  L (k,a)Z (k)]P(Z,t) (4.66)
Oa ij i at a aa 8
62
and the diffusion term is
J1 (Z,t) = R (k) () P(Z,t). (4.67) The matrix R a(k), the diffusion matrix, is just the correlation matrix of the stochastic forces. Equations (4.65)  (4.67) are the usual results for the relationship between a linear Langevin equation and the corresponding FokkerPlanck equation.
CHAPTER 5
NONLINEAR FLUCTUATING HYDRODYNAMICS
In this chapter, we will consider the consequences of the
nonlinear terms occurring in the general fluctuating hydrodynamics model. We begin by considering a naive form of renormalization which illustrates in an intuitive way the concept of renormalization. Using this method, we will derive the lowest order renormalization of the shear viscosity explicitly deriving the nonanalytic dependence on shear rate. In the following sections, we will present a more formal development of renormalization which will result in a selfconsistent formulation of the nonlinear problem. The chapter concludes with the connection between the naive and formal renormalization.
Naive Renormalization
Consider a simple, one dimensional, nonlinear, stochastic equation of the form
a 3
Sx + LOx + N x f(t) (5.1)
with a Gaussian, white stochastic source f(t) having correlations,
f6 (t t2) (5.2)
63
64
and all higher order cumulants equal to zero. Furthermore, let the average of x(t), , be denoted by x0(t). Then, the equation of motion for x0(t) is found by averaging eq. (5.1),
a Xo(t) + Lox0(t) + N = 0. (5.3)
If we write x(t) = x (t) + 6x(t), then (5.3) becomes
a X0(t) + [LO + 3NO<(6x(t)) 2>]x0(t) + No0x(t) =  No<(6x(t)) >.
(5.4)
We see that the equation for xo(t) has the same form as that for the fluctuating variable, x(t), except that the "bare" coefficient Lo has been renormalized by the nonlinearity to become the "dressed" or renormalized coefficient L,
L = L0 + 3N0<(6x(t))2>. (5.5)
If we subtract eq. (5.4) from (5.1) and solve the resulting equation for dx neglecting the nonlinear term, and use this to compute L, we find that
L = LO + 3N f /2L + o(N ) (5.6)
where we have imposed stationarity,  = O, to fix the atequal time correlation functions.
equal time correlation functions.
65
The previous calculation can be reformulated in terms of
correlation functions. The equation of motion for the twopoint function, <6x(t)6x(0)> s C(t), is found from (5.1) to be
at C(t) + LoC(t) + N0<(6x(t))36x(0)> = 0 (5.7)
where we have used the condition  0. To lowest, NO = 0, order, we find that the Laplace transform of C(t), C(z), is
C(z) = C(0)/(z + LO) + o(N ) (5.8)
which has a pole at the bare transport coefficient L0. If we calculate the correlation function in (5.7) to lowest order, we find that
t ddd3f0 dT4 0(tT 1) L0(t2) L (tT3) L0 4 <(6x(t)) 6x(0)> = e e
+ o(N )
t 0 e2L (t ) L (tT2 ) = 3 jtd j d'r2 e e f + o(No)
* 1  2 0 0
1 1 L0t 2
=3 2L0 2L0 e fo + o(No)
0 (0)
C (t) + o(N0). (5.9) 2L0 0
66
So,
SC(t) + (L + 0) C(t) = o(No2) (5.10)
t 0 2 L0 0 and
3 2
C(z) = C(O)/(z + LO + N f /LO) + o(N ) (5.11) which has a pole at the "renormalized" coefficient defined in (5.5). This coincidence is easy to understand. Suppose that at time t=O, the distribution function, p(x), is perturbed from the steady state function, ps, in the following way
h6x
p(t=0) = ps e . (5.12) The auxillary field, h, might be chosen such that h * ss = XO i.e., to yield some particular initial value for the macroscopic variable. For small deviations from the steady state (i.e., h small so h  ss is small), we find that
 x0 = C(t)h + o(h2)
 (C(t)/C())[h  x 0] + o(h2). (5.13) Thus, the correlation function acts as the Green's function for small
67
deviations of the macroscopic variables. They must, therefore, share the same linear relaxation times. A general .proof of this property for Markovian processes has been given by Dufty et al.45
In the same way, we can derive an expression for the renormalized shear viscosity. If we take the average of eqs. (2.45), we require that we obtain the macroscopic hydrodynamic equations. For shear flow, the only nontrivial equation is the heat equation,
8 2
e  nRa = 0 (5.14)
where nR is the renormalized viscosity. Comparing (5.14) with the average of eqs. (2.45) yields
2
n = a
R ii Dr 1i r tij
2
= a
2 8
= a < ui(r,t)] p(r,t) u (r,t) u (r,t)> (5.15)
where we have used temporal and spatial translational invariance to set overall derivatives of the correlation functions to zero. The third line of (5.15) follows directly from the momentum equation.
Using our results from the previous chapter, eq. (5.15) can be written, to lowest order in the nonlinearity, in terms of the equal
68
time correlation functions,
nR = a p 0
(2w) 6d k[k k y 33(k;a) +k k y44 (k;a)  k 45(k;a)]. (5.16)
This equation for the lowest order renormalization of the viscosity is an example of mode coupling. That is, the original nonlinearities in the equations of motion have coupled two different hydrodynamic modes, analogous to phonons in solid state physics, and the correction to the viscosity is due to the interaction of these modes, analogous to the phononphonon scattering contribution to the coefficient of thermal diffusion in solids. The coupling of the modes in (5.16) is buried in the equal time correlation functions. It can be seen by referring to eq. (4.38). Contributions also exist due to three mode coupling, four mode coupling, etc.
Upon evaluation, we find that
1/2
nR n0 + nl/2a (5.17) with
S=  kBT 2.56 + 4.(5.18) (2vB)3/2 (r0)3/2
Thus, it is seen that a nonanalytic dependence of the viscosity on shear rate arises as a result of mode coupling. This is in contrast to the analytic dependence one would obtain from the standard ChapmanEnskog solution of the Boltzmann equation (or the Enskog equation; see
69
Chapter 6). This result agrees with that of Ernst et al. and with that of Yamada and Kawasaki (although their numerical evaluations are incorrect). Thus, this result indicates the equivalence, with the present assumptions, of this model and the others appearing in the literature.
Formal Renormalization
In this section, we wish to develop a more general method of
renormalization than the naive scheme, allowing the calculation of all renormalized properties. To this end, we begin by formulating the nonlinear model in a more general notation,
 Y + vmEy] + D [y]F [y] = M [y]n (5.19)
where v [z] represents the nonlinear Euler terms in eq. (3.23) (i.e.,
m
the terms V * j = V * (pu) in (3.23) and eva and v j8 in (3.24)). The matrix Dmn[Y] is called the diffusion matrix and F ny] are the thermodynamic forces which are the variables thermodynamically conjugate to the ym's. The product DmnFn represents the dissipative terms in (3.24). The reason for returning this form of the equations is that in their derivation of fluctuating hydrodynamics, Zubarev and Morozov5 show that thermodynamic consistency imposes the fluctuationdissipation theorem,
D mn[y] = M m,] Mnn[y] Amn,; (5.20)
70
"thermodynamic consistency" meaning that the equilibrium distribution function is a solution to (5.17).
Equation (5.17) is highly nonlinear making a general analysis
difficult. However, in equilibrium it has been found that truncating eq. (5.17) at second order in the nonlinearities provides a qualitatively correct description of such nonlinear phenomena as long time tails (i.e., the algebraic decay of equilibrium correlation functions). In fact, in equilibrium, the dominant mode coupling effects arise from the Euler terms in (5.19). We shall, therefore, consider only these nonlinearities in the rest of the discussion. It can be shown that, to order (ak), this is equivalent to taking both Dmn and Mmn to lowest order (i.e., Dmn[y]  Dmn[Y 0]) thus preserving the fluctuationdissipation theorem, eq. (5.20), in this truncated model of eq. (5.19). Furthermore, Fn[y] is expanded to first order in the fluctuations in order to preserve the usual linearized transport laws.
With these approximations, and expanding the variables ym about their steady state averages y0m, eq. (5.19) takes the form,
z (k,t) + L (k,t) Zb(k,t) + eV (k,k k )(z (kl t) z (k ,t) t a b ab c 1 2 b 1 c 2
 )
= Ra(k,t) (5.21)
where L = a k  + H (k;a) is the linear operator studied in
ab iji k ab
71
the previous chapter and za  Ya  YOa (details are given in Appendix B). The vertices, which arise solely from the Euler terms, have the form
Vab (k,k ,k ) = A (k,k ) 6(kk k ) (5.22)
abc 2 abc1 1 2
where the explicit form of Aabc(k,q) is given in Appendix B. The constant E is introduced as a notational convenience in the development of a perturbative expansion in the nonlinearity.
Equation (5.21) can be rewritten using the Green's function introduced in Chapter 4 as
(0) r
z (k;t) = z( (k;t) + eft d{G (kql ;tT) Vabc(q,,3) a a aa 1a'bc 1 293
x (zb(q2,T) zc(q3,T)  )} (5.23)
where a sum over repeated wavevectors is implicit and where
za(0)(k;t) Gaa,(k,q;t) za,(o,0) + ft d Gaa (k,q;tT)Ra,(q,T).
(5.24)
Using (5.23), a perturbative solution to eq. (5.21) is easily constructed.
From eq. (5.21), we can immediately construct the equations of motion of the twopoint correlation functions,
72
Ca (k,k';t) = ,
C ab(k,k',t) + L (k;a) C,(k,k';t) + M (k,k';t) = 0 (5.25)
St ab aa a'b ab where we have defined
Mab(k,k';t) = Vacd N,1,q2) (5.26) and have used the fact that = 0 for t'
Mab(k,k';t) ftdT aa (,k ;t) Cab(k ,k';r). (5.27)
a,0 aa 1 alb 1
The significance of this is that eq. (5.24) becomes a linear equation for the correlation functions and effective viscosities, etc. as would be measured in the laboratory may be identified.
Laplace transforming (5.24) in conjunction with (5.27), one obtains
z Cab(k,k';w) + Laa(k;a) k,k';w) + L ,(k,k1;w) Cab (k ,k';w) ab aa, alb ' w aa' '1 alb ;'1 w
= Cab(k,k';t=0) (5.28) and the Fourier transformed selfenergy may be computed from (5.27),
73
ab("k' ;w) M(k,;w) C (,cb( j';w). (5.29)
As discussed in the previous chapter, Cab(k,k';t) = Cab(k,t)6(kk'(t)) and it is easily verified that, order by order in perturbation theory, Mab(k,k';t) m Mab(k,t)6(kk'(t)) which requires that lab(k,k',t) m Iab(k,t)6(kk'(t)). Thus, eq. (5.28) becomes
Cab(k,t) + D(k)C (k,t) + H (k;a)C (k,t) + eM (k,t) = 0 t ab ab aa a'b ab (5.30)
where we have introduced D m a jk.i . The equation for the selfenergy becomes
Ma (k,t) = fd aa (k;tT)Ca, (k(Tt),t))
rt D(tr)
= dr (k;tT)e C (k,T) (5.31)
0 aa' a'b
which shows the selfenergy to be an operator in kspace. Laplace transforming (5.31) yields an expression for the selfenergy,
wt Dt  1 fdt e ,(k;t)e = M (k;w)Cb (k;w).
aat (k,w) (5.32)
74
and the Laplace transform of (5.30) is
Sd ; )e (w+D)(T) (
[w+DICab(k,w) + H (k,a)C ab(k,w) + efo taa (k;())e C (k,w
 Cab(k,t=O). (5.33)
Equations (5.32) and (5.33) are the primary results of this
section. The linear mode coupling which was the unique feature of the linear problem is seen to play an important role here as well. The operator nature of the selfenergy is a result of the linear mode coupling and expresses the fact that new gradient couplings, as well as the expected new hydrodynamic couplings, arise as a result of the process of renormalization. This complication is, however, removed in the context of perturbation theory in e because we have an expression, to lowest order in E, for the operator (w+D) acting on the correlation functions (c.f. eq. (5.33)). It is seen that the selfenergy function shifts the poles of the correlation function renormalizing the transport coefficients as in eq. (5.11).
Application of Formal Renormalization
We now wish to apply the renormalization formalism to the
calculation of the renormalized shear viscosity and, in so doing, to make a connection with the naive renormalization method. As a simplification, we will restrict our attention to the contribution arising from the shear modes. In this case, the irreversible part of the stress tensor, tij, can, in general, be written as
75
t . j(5.34) ij ijm Jrl 5
where Vijlm is the viscosity tensor which, in the unrenormalized theory, has the NavierStoke's form,
S = v(6 6 + 6 6 2 6 6) (5.35)
ijlm il jm im6jl 3 1 6ij1m) ij 1m
where v and K are the usual kinematic and bulk viscosities, respectively. The macroscopic equation for the energy density is
 1 1
a e + r [(e + p)p j + t jl + s] = 0 (5.36)
For macroscopic shear flow, using eq. (5.34), the heating equation becomes
Se + ij aam 0. (5.37)
 e+ jlm ijalm
Thus, the renormalized viscosity calculated earlier from the heating equation corresponds to the yxxy element of the general viscosity tensor.
We now calculate the selfenergy to second order in the
nonlinearity (i.e., to second order in E). Using eqs. (5.23) and (5.26), we find that
76
Ma (k,k';t) EVad(_,q,2)
+ 2E 2dVacd (k,( l 2)Gce(a, '3 tT)Vefg(434'9 5)
x <(0) (0) (0) 3
x } + o(3 ) (5.38)
where the factor of two arises from the symmetry of the vertices. Using eq. (5.24) and causality, this reduces to Mab(k,k';t) = EVacd(k,ql92.2)Gcc,(1',3;t)Gdd,(2'q4't) x 0+ 2E2od~Vacd , 1,2)Gce 1'9 3 tT)Vefg 3'44'95) G f,(4 ,q4;T)Ggg' (5'9;T)Gdd'( 2';t) x } + o(3).
(5.39)
In equilibrium, the equal time fluctuations are Gaussian distributed making the o(c) term in eq. (5.39) zero and allowing easy evaluation of the o(E ) term. In the steadystate, the equal time fluctuations are not necessarily Gaussian distributed. Stationarity, in the steady state, allows us to write
77
za(k;O)  z)(k;O) + ef0 d Gaa ' (kq'T )Vabq1' .2'3) x (zb(q2;T)Zc(3;T)  )} (5.40)
where
)(k;0) = d G (k;)R ) (5.41)
a   aa' ; )Ra 1T) and
(0) (k,t) = td G aa(k,q;tT)R (q1r) (5.42) where in (5.42), t>0. Thus, the equal time fluctations may be determined perturbatively and are Gaussian distributed to zeroth order
in :.
Equations (5.39)  (5.42) provide us with enough information to calculate the selfenergy function from eq. (5.28). The result is that
aa ,(k,k";w) = e cdt ewt IV (k,ql )cc ( 1'3;t)Gdd' (2'4;t) x ba (k',k";z)
78
4+ 4c2dt e  Vad(t ,l1,2)Vega,(3 , k4,k")
x Gde (2,3;t)God21,9;t)
+ o(E3). (5.43)
To recover our previously calculated lowest order mode coupling contribution to the viscosity, we first replace the equaltime correlations in (5.43) with their local equilibrium values. The term proportional to c in (5.43) then vanishes and (5.43) becomes
 kBT
E (k,k";w) = fdt e A ,(k(t)q(t),q(t)) aa' PO acd eya
D(k)t
x Gde(kg;t)Gcd (;t)e  }D6(k+k").
(5.44)
The momentummomentum selfenergy upon evaluation is found to have the form
2 (k,z)  kik j(k,z) (5.45)
which identifies it as the renormalization of the viscosity tensor.
(0)
If we write ijl = Vijlm + 6v'ijlm then the long wavelength, long
jlm viscositym appearing in eq. (5.37) ism
time hydrodynamic viscosity appearing in eq. (5.37) is
79
6v = lim (k,w) (5.46) yxxy k+O yXPxPy
w+O
from which the "naive" result is recovered.
The replacement of the static correlation functions by their local equilibrium values requires some comment. For a stationary process, the correlation of the stochastic sources is related through fluctuationdissipation theorems to the static correlations. The mode coupling contributions to the static correlations may thus be interpreted as a renormalization of the amplitudes of the stochastic sources. By neglecting the equal time threepoint function in (5.42), we were left with the propagator renormalization.
There is another way to understand the replacement of the static correlations in the present calculation by their local equilibrium values in order to find agreement with the "naive" calculation. The only information about the "statics" that was used in the "naive" calculation was the stochastic source correlations. But, as explained in Chapter 3, these are approximated, in the unrenormalized theory, by their local equilibrium values. It is, therefore, reasonable that the same information must be used in both calculations to obtain agreement. We speculate that if the first order renormalized forces were used in the naive calculation, one would obtain the same result as is obtained here without the local equilibrium replacement.
CHAPTER 6
SHORT TIME MODEL
In this chapter, we discuss the explicit form of the short time model of hydrodynamic fluctuations introduced in Chapter 3. We recall that this model is intended to describe hydrodynamic fluctuations on length scales as small as atomic lengths. We begin by discussing general properties of the distribution functions. This is followed by a derivation of explicit expressions for the necessary matrix elements in terms of the reduced distribution functions. The final two sections describe the approximations introduced to evaluate these matrix elements and the explicit form of the resulting equations governing the time evolution of hydrodynamic fluctuations.
General Properties of the Reduced Distribution Functions
In this section, we introduce two properties of the distribution functions which will be of use in determining the general nature of the expressions, eqs. (3.31)(3.34), to be evaluated below and which will aid in their evaluation. These properties are the stationarity, in time, and the modified translational invariance of the distribution functions.
As discussed in Chapter 2, an external nonconservative force is applied in computer simulations of sheared fluids to counterbalance viscous heating of the system. The significance of this is that it
80
81
yields a steady state. If we denote by pss(r) the distribution function (density matrix) of the steady state system, depending on all of the phase space variables r i then we can write the steady state Liouville equation as
0 a p (r) + ï¿½p(r) (6.1)
= ï¿½ p(r) (6.2)
where ï¿½ is the Liouville operator. The explicit form of ï¿½ is generally given by
A N Vk 1 a
= 1 V( nn' (6.3)
n m nï¿½n "n 2 "n an ext
n=1n nn' n 2n
where V( nn ,) is the pair potential. For hard spheres, the potential is singular rendering eq. (6.3) ill defined. As discussed in Chapter 3, however, eq. (6.3) can be replaced by a well defined expression which produces a statistically equivalent dynamical system
^ N ^
ï¿½ = n + 1 [T(n,n') + T_(n',n)] + ï¿½ (6.4)
n=1 n nsn'
where we have included a term representing the externally applied thermostat alluded to above. The explicit form of the T operator and its various adjoints and timereversed forms are given in the literature and will not be repeated here. The point we wish to make
82
is that (6.2) allows the ï¿½+ operator to be moved in averages, i.e., for any two phase functions, A(r) and B(r), it can be shown that
= <(ï¿½ A(r)) B(r)> . (6.5) s  ss
The subscript on fï¿½ refers to the fact that T+ goes to its time reversed form T..
Another important property of the distribution functions is their modified translational invariance. Clearly, because the local flow velocity varies spatially, the distribution function cannot be translationally invariant. That is, if
i *
pi,aqi , ,'i  R (6.6)
for an arbitrary displacement R , then
Pss(r) * Pss(r*). (6.7)
However, if we make the transformation
i' 4 , pi  maij Rj i  R , (6.8)
then we expect the distribution function to be invariant since the macroscopic state is characterized by the local velocity and the local velocity gradient and, while (6.6) does not preserve these properties, (6.8) does. Mathematically, it is easy to verify that (6.8) is
83
consistent with the LeesEdwards boundary conditions while (6.6) is not. We, therefore, look for solutions of Liouville's equation which are functions of the variables
Pj s pi  maijqj
q q  q (6.9)
which are invariant under (6.8), i'* i i R'2j. (6.10) Henceforth, we will write Pss(r) P, p({p, qij).
Explicit Form of Matrix Elements
The short time model is given in Chapter 3 as
 y (k,t) + L (k,k';a)6y(k',t)  0 (6.11) where
L (k,k';a) = <(ï¿½y a(k))6yy(k")>ss (g )y(6.12) and
84
g8(kt",k') = <6y (k")6y (k')> ss (6.13) The variables, 6y (k), are again taken to be the fluctuations of the conserved densities around the steady state:
^ ^ 1 ik*
6y1 (k) = p(k) = 1 me iqn
6y2(k) = 6e'(k) =  k ek
2 2m 2 Be
( = 6j'(k)  . ( eikp)n (6.14) V n
which are summarized as
6ya(k)  Ya (') e n, (6.15) /V n
where V is the volume of the system.
We now define the reduced distribution functions,
nN! fdX' ... dX' p (r') (6.16) nnf(n)X''X~) (Nn)! n+1 n ss
where X' (n,2n) and we note that the Jacobian of the change of variables X + X' is unity. Using (6.16) and the independence of Pss on the center of mass, we find that
ga (k,k') = (2)3 V 16(k+k') g_ (k), (6.17)
85
where
ga8(k) = n fdp f() 1 a Y8
1 2 (2) ik*r
+ 2 n d dpdr2f (1 ,, )12Yy 1 ()y ()e 12.
(6.18)
The explicit form of ï¿½ +y(k) is easily found to be
^ ^ 1 k En Y Eik*n
ï¿½ ya(k) = i Zk  Y, ()e a a n
v n n
+ .  [T+(n,n') + T+(n',n)] Y(n)ei 'Q
/V nln'
 k a j y (k) + ï¿½ y (k). (6.19)
i ij3ak a ext aFrom this, it follows that
<(ï¿½+y (k))6y (k')> = (21r) 6(k+k')M (k)
 k.a ~ Sij k a 8 ss + <(ï¿½e y (k))y (k')> , (6.20) where
86
M (k) = ink.fdp 1 Ya(()Y1 ()f ()
2 (2)
+i kd2 m a 1 Y )Y )f(2 ;12)
1 (2)
+ n2f dpdd 2( T+(1,2)+T+(2,1)])Y (aY 8 (2{' 12)
2 ik'r
n  12
d+ 2fd dd 1 2([T +(1,2)+T (2,1)])Y 1 (n )Y (2)e
x f(2)
2 ^ ik*r
+ fdpdpdp dql2dq_3(E +(1,2)+T (2,1)])Y (p1 ))Y(3)e
x f(3) (1' 2Rj~I;12'13). (6.21)
Approximations
The evaluation of eqs. (6.21) and (6.18) require knowledge of the one and twobody reduced distribution functions for the steady state. In general, these functions cannot be determined explicitly and, so, we introduce two important approximations.
The first of these is the neglect of momentum correlations, which is to say we assume the twobody distribution to have the form
f(2) 2 (2) 2
Sf() ()f )g(12), (6.22)
87
where f(1)
where f (p) is the onebody reduced distribution function, which cannot depend on 1 due to the modified translational invariance, and the function g(q12) is the steady state analogue of the pair correlation function (i.e., g(q12)dq12 is the probability of finding a particle, 2, in the volume da12 if there is a particle, 1, at position .11). In equilibrium, this approximate form becomes exact. However, away from equilibrium, we expect that particle velocities are indeed correlated (i.e., (6.22) implies that ss = 0 but in general this need not be the case). This approximation, which is analogous to Boltzmann's Stosszahlansatz,43 is standard in the development of kinetic equations. In particular, it is used to derive the Enskog equation which is known to give a qualitatively accurate description of fluids near equilibrium.43 With this approximation, the threebody term in (6.21) only contributes for 8=1. In Appendix D, it is shown that this term can be written, using an exact identity following from stationarity, in terms of twobody functions.
Having adopted the form (6.22), we must still determine the onebody distribution and the pair correlation function. The second approximation we use is to replace the steady state pair correlation function, g(a12), by the local equilibrium pair correlation function, 8LE(I 121) ' gLE.(q12). The pair correlation function is related to the densitydensity equal time correlation function and knowledge of one of these allows determination of the other. This approximation, which, like the previous one, is used in the derivation of various kinetic equations, may be thought of as a first guess to be "renormalized" using the density autocorrelation function calculated
88
from the resulting equations. This approximation then becomes the first iteration of a selfconsistent calculation. In fact, these two approximations require that f(l)( P) satisfy an Enskog equation which is the usual dense fluid analogue of the Boltzmann equation.
With these approximations, our model is completely defined. The onebody distribution can be obtained from the first equation of the BBGKY hierarchy (i.e., by integrating eq. (6.2) over coordinates X)...XN and using eqs. (6.4) and (6.22)):
m (a f ( 1') = fdX2T (12)f )f(f () (2)gLE12)
+ sources. (6.23)
The sources in eq. (6.23) arise from the nonconservative forces in the Liouville equation (eq. 6.4). However, it can be shown that they contribute to second order in the shear rate and so are neglected in the present calculation which is only carried out to first order in the shear. While the previous approximations served to construct a model, we now make an approximation with respect to solving that model; we will only work to first order in the shear rate. Thus we expand f(,
f ( ) = fLE , )(1 + a*( ) + ...) (6.24)
where fLE(pl) is the local equilibrium distribution function for uniform shear flow,
89
fLE(p) = n(2mnkBTo )3/2 exp  (pl /2mkBTo). (6.25)
With these approximations, one is able to determine ï¿½1(P) in eq. (6.25). As this result exists in the literature, we only quote it here:
a( (P) = ajnk I/rO(P(iPj  1 6 )
1 4
nk = [1 + j5 n*X]nB
nB = 2 (m//80) (6.26) 16a
and X m gLE( ). The constant, nB, is the kinetic part of the viscosity as calculated from the Boltzmann equation. It should be noted that we have only sought the steadystate solution of eq. (6.22).
With these approximations, the matrix g a(k) may be determined to be
g11(k,a) = nS(k)
g22(k,a)  n(k T )
ga (k,a) = PokBT 0 a  mnk(a. +a. ) a=i+2, 8=j+2=3,4,5
gaB = 0 all others (6.27)
90
while the matrix M a(k,a) can be written as
M =n<  (') YB(E')>
+ 6a2 n aij((P j 3 p'261j)/m, y 8('))
+ Ia8  Sa8 (6.28) where <"'>E denotes an average over fE( ') truncated to first order in a. The collisional contribution, l.8, is
I g(k,a) = n Xd'dRdp l2 Y(e) T(12)f (plf (2) ^ ik*l12 ik12 (6.29) x [Y B(g) + Y B()e  1e (6.29) and the source term, S8a, is
SaB = <(ext (k))[k)) )  k  ss]> ss. (6.30)
2
It is shown in Appendix D that S  a, allowing us to neglect it here.
Hydrodynamic Equations
Once the various integrals have been performed, we find the
following form for the short time hydrodynamic equations, where the
91
stars indicate dimensionless quantities with lengths scales to a and times to tE,
* 2
DtP  ixk* = 0
* * * 2 * ^^ * Dt6e + [D T(X)x + a c k k ]6e
*
^ * ^ ^ * ^ ^ ^ *
 i[xh(x)k  a (c + y)(k 6 + k 6 )  a k k k a]6j
a 2 y xy 3xya a
X 2
= o(a2)
Sa a*x ^ ^ D t6j a  i[xP (x)k  a (k 6 + k 6 )]6p t a p a 6yS(x) y ax x ay
i[xP (x)k  a (c4 x)(k 6 ax+ k 6 )]6e
e a +y ex x sy
2 * * 2 *
+ v (x)x 6 + (v' (x)v (x))x k k ]6j
+ a [c6(6ax 6y + 6 ay6 ) + C kx k y6
+ o8k (k6x + k 6 ) + (k ky6 x + kx6 ay)k8
^ . ^ * * * 2
+ I0 kxkyk k $16j + a 6 6 6j = o(a ) (6.31)
where x = ko and,
92
* 2 2 DT (x) = (1  Jo(X))/x
2 2 v' (x) = (1  J(x) + 2j2(x))/x2
2 2
v (x) = (1  j0(x)  j2(X))/x
* It
h (x) (1 + 3yj1(x)/x)
36y2
* T S1x)
P (x) S (x)
36y
P (x) (1 + 3yj (x)/x) (6.32)
e 3 y1(x)3yj
are wavevector dependent generalizations of the thermal conductivity, viscosity, shear viscosity, enthalpy, density derivative of pressure and energy derivative of pressure. The constants y = 2 n
* 5w 2 *
and v [1 +  y]; v is a dimensionless form of the kinetic
24y 5
contribution to the viscosity, nk. The functions jn(x) are nthorder spherical Bessel functions with j1(x) = (sin x)/x. The functions c = ci (x;a) are given in Table 6.1. Finally,
* 3 .* 2 * a * a p pa j = jt E/a, etc. while Dt aijki . is the at j
convective derivative in terms of dimensionless time t* = t/tE and a = atE.
At this point, we make a few comments about the general nature of these equations. The first of eqs. (6.31) is the familiar continuity equation describing the conservation of mass. The second two
93
equations describe, respectively, the conservation of energy and momentum. However, due to the fact that the effects of atomic structure are included, they are more complicated than at small wavevector. These complications disappear if we take x << 1 and expand to second order in x in which case we get a short time version of the linearized NavierStoke's equations. We note that in h, Pe' Pp and c2 we have picked up second order contributions to the pressure in the long wavelength (small x) limit. To illustrate the meaning of the new couplings, we have, in Table 6.2, given the expansions of the coefficients ci to second order in x. It is seen that all of the collisional contributions to the nonequilibrium couplings are of third order or higher in the gradients except for c2 (recall that these coefficients are all multiplied by the shear rate, a, which is a first order gradient). Thus, they do not appear in the ordinary NavierStokes equations but rather in the higher order, so called Burnett, equations. They are always present in the analysis to be presented below, however, as we keep all orders in ka.
Also, these equations contain wavevector dependent "transport
coefficients," eq. (6.32), which are the high frequency limit of the general wavevector and frequency dependent transport coefficients,
DT(X) = lim DT(x,w), (6.33) W).
etc. Wavevector dependent susceptibilities such as the pressure derivatives also appear. An important feature of the transport
