A numerical test of stress correlations in fluctuating hydrodynamics
AA numerical test of stress correlations in fluctuating hydrodynamics
Michael Schindler
Laboratoire PCT, UMR “Gulliver” CNRS-ESPCI 7083, 10 rue Vauquelin, 75231 Paris Cedex 05
Abstract
The correlations of the fluctuating stress tensor are calculated in an equilibrium molecular-dynamics simulation of aLennard–Jones liquid. We define a coarse-grained local stress tensor which can be calculated numerically and whichallows for the first time to determine the stress correlation function both in time and in space. Our findings corroboratethe assumptions made in fluctuating hydrodynamics as long as the liquid is isotropic, that is in bulk. In the vicinity of arigid plate, however, the isotropy is restricted, and major modifications must be done with respect to the usual theory.Among these are the appearance of five different viscosities instead of two and a non-trivial dependence of the distancefrom the wall. We determine these viscosities from the simulation data and find that their values are very different fromthe bulk values. We further find much longer relaxation times of the stress correlations than in bulk.
1. Introduction
The theory of fluctuating hydrodynamics starts withclassical hydrodynamics, which describes the average mo-tion of a fluid and which is formulated in terms of smoothtensor fields, and enriches it by additive fluctuations. Itcan be understood as the attempt to bridge the gap be-tween microfluidics and nanofluidics from the large-scalesside. The central question is of course, what propertiesthe fluctuations must have in order to correctly representthe microscopic states of the liquid in a given problem?The standard approach [1, 2] is to split the full veloc-ity and stress tensor into a macroscopic (denoted by anoverbar in the following) and into fluctuating parts (de-noted by a tilde). The linearised governing equations takethen the form of a two-scale process, with the standardcompressible Navier–Stokes equations for the macroscopicvelocity ¯ v ( x , t ) and the corresponding Newtonian stresstensor ¯ σ ij ( x , t ),¯ σ ij = − δ ij ¯ p ( x , t ) + η ( ∂ i ¯ v j + ∂ j ¯ v i ) + (cid:16) λ − η (cid:17) δ ij div ¯ v . (1)Here and in the following, η denotes the shear viscosity ofthe fluid and λ its volume viscosity. As we will be con-cerned only about equilibrium situations, the macroscopicequations reduce to¯ ρ ( x , t ) = ρ , ¯ p ( x , t ) = p , ¯ v ( x , t ) = , (2)with constants ρ and p . The fluctuating fields satisfy alinearised Navier–Stokes equation with additional randomnoise s ij , ρ ∂ t ˜ v = div(˜ σ + s ) , (3) ∂ t ˜ ρ = − ρ div ˜ v . (4) Email address: [email protected] (MichaelSchindler)
The stress tensor in equation (3) is again split into twoparts, the first of which satisfies the same Newtonian def-inition (1) as the macroscopic stress tensor, only in termsof a fluctuating velocity ˜ v ( x , t ) and of a fluctuating pres-sure ˜ p ( x , t ). The second contribution is an uncorrelatedfluctuating stress tensor s ij ( x , t ). The tensor ˜ σ ij is meantas a coarse-grained stress which is fluctuating, but not asrandomly as s ij , see for example the book by Kubo et al.[3] for coarse-graining concepts, and Ref. [4] for substi-tute Markov processes. The same two-step separation hasbeen used by Hauge and Martin-L¨of [5], Chow and Her-mans [6], but without commenting on why it is necessary.Equation (3) has the structure of a Langevin equation,with s ij the random driving. This random stress tensorvanishes on average, (cid:104) s ij ( x , t ) (cid:105) = 0, and it has to satisfythe fluctuations–dissipation theorem of the second kind . Itis thus uncorrelated in space and in time, (cid:10) s ij ( x , t ) s kl ( y , s ) (cid:11) = 2 kT A ijkl δ ( x − y ) δ ( t − s ) . (5)Here, (cid:104)·(cid:105) denotes the average over different realisations, ora time average; k is Boltzmann’s constant, and T the tem-perature of the system. The correlator is further assumedto be isotropic, which fixes the fourth-rank tensor A ijkl to [2] A ijkl = η (cid:2) δ ik δ jl + δ il δ jk (cid:3) + (cid:16) λ − η (cid:17) δ ij δ kl . (6)This special form of A ijkl follows from the fact that itmust be isotropic, thus invariant under any rotation of thecoordinate system. It is therefore a linear combination ofthe only three isotropic tensors of rank four, which are δ ij δ kl , δ ik δ jl , and δ il δ jk [7]. The factors of proportion-ality are found to be the three viscosities of a fluid [8].The additional symmetry of the stress tensor, leading to A ijkl = A jikl = A ijlk eliminates the rotational viscosity. Preprint submitted to Chemical Physics October 24, 2018 a r X i v : . [ c ond - m a t . s o f t ] M a y n the vicinity of a rigid boundary, however, isotropy isnot granted in the same way as in the bulk. Instead, thestress correlator should reflect the fact that the normaldirection is different from the two tangential directions.In this paper, we therefore ask which tensorial form thecorrelator will adopt when evaluated close to a rigid wall.In order to start with the most simple situation, we regardthe stress correlator in thermodynamic equilibrium.We must further ask, what “close” means, that is upto which distance from the wall the different symmetry issensed by the fluctuating part of the stress tensor. Nearthe boundary, we thus question also the spatial delta func-tion in Eq. (5). The problem with the spatial delta func-tion becomes apparent when applying fluctuating hydro-dynamics for example to the Brownian motion of a spher-ical particle. This has been done in the literature on theautocorrelation of a Brownian particle [5, 6]. If we naivelyintegrate the random stress tensor over the particle in or-der to obtain the random force, and then construct theautocorrelation of that force, the result would always beinfinite. The reason is that we then integrate the three-dimensional delta function from Eq. (5) only over the two-dimensional surface of the particle. The third dimensionthus remains a delta function, evaluated at zero.The spatial delta function must be questioned also froma less formal point of view: The standard argument for theabsence of spatial correlations is based on thermodynam-ics [1]: The fluid volume is cut into several small but finiteportions ∆ V . In each of these volumes, the entropy pro-duction is calculated as the volume-integral of the tensorialreduction of the stress with the rate-of-strain. The fluctua-tions of the stress tensor is then taken from the variationsof the entropy production around the equilibrium state.As this integral does not contain any surface contribution,the fluctuations of two stress tensors attributed to two suchboxes are independent of each other. This property is thenkept even in the limit of infinitesimal box size . It appearsas if there should be a lower spatial bound to the validityof this argument, namely when the volume size become ofthe order of the distance between the molecules of the liq-uid. As thermodynamic arguments are used, the thermo-dynamic limit of many interacting molecules is assumed tobe done before the box size can be sent to zero. It appearsthus as a reasonable question to ask, down to which spatialscale does the delta function yield a good description? Itwill turn out in the following, that it works amazingly well:In bulk fluid, the delta function holds down to distancessmaller than the intermolecular spacing.It is the aim of this paper to test the validity of Eqs. (5)and (6), describing stress correlations—in particular, tofind deviations in the vicinity of a rigid surface. In or-der to do so, we start from a fully microscopic simula-tion of a liquid and extract the stress tensor using coarse-graining. As the simulation is done at thermal equilib-rium, the macroscopic fluid velocity ¯ v vanishes, leavingonly the thermodynamic pressure in the macroscopic stresstensor ¯ σ ij = − p δ ij . We correlate the coarse-grained stress tensor at different positions and at different times in or-der to verify both the isotropy in Eq. (6) and the delta-functions in Eq. (5). This will be done both near a rigidwall and in bulk fluid.The precise form of the stress tensor correlations is po-tentially important for all fluidic systems including bound-aries. In microfluidic applications, and even more in nano-fluidic ones, the role of the boundaries is essential [9, 10,11, 12]. In nanofluidics, the presence of boundaries evenintroduces different time scales in the relaxation dynam-ics [13]. For example, Brownian motion of a particle im-mersed in a viscous fluid, which is one of the most classi-cal systems both in the hydrodynamic and in the stochas-tic community, is in most cases a microfluidic problem—depending on the size of the particle. Both for large andfor small particles, the velocity fluctuations are known toexhibit long-time tails in their autocorrelation-function.The theoretical explanations of this phenomenon, as faras they consider both the particle and the fluid, assumethe isotropic form of Eqs. (5) and (6) for the fluctuatingpart of the stress tensor—even when this stress tensor isintegrated over the surface of the particle to yield the totalfluctuating force acting on it [5, 6]. The investigation inthe present paper is one of two necessary steps to renderthe theoretical description of the long-time tails more con-sistent. We will find below that it is not the usual viscositywhich occurs in the formula of the force autocorrelation,but a viscosity which is a property of the fluid and thestructuring effect of the rigid particle itself. The otherstep will be published in Ref. [14].We will further find that the inclusion of boundariesmakes it necessary to pass from homogeneous viscosities(which are defined in the limit of homogeneous space with-out any scale) to rather local viscosities. Such viscositiesreflect more details of the geometrical conditions and ofthe time scales involved.
2. Symmetry considerations
The form (6) of the stress correlator follows directlyfrom symmetry considerations. We will now apply thesame considerations to derive the corresponding expressionfor a liquid close to a wall. Similar formulae are foundin the literature on nematic liquids [15]. The followingsymmetries apply to the stress correlator:(A) The most rigorously enforced symmetry is that the stress tensor is symmetric . This symmetry is inheritedfrom the microscopic definition of the stress tensor and isthus not subject to random fluctuations. It leads to theconditions (cid:104) s ij s kl (cid:105) = (cid:104) s ji s kl (cid:105) = (cid:104) s ij s lk (cid:105) . (7)By this symmetry, the 81 components in the full corre-lation tensor of rank four in three-dimensional space arereduced to 36 different ones.2 + 43 η η η + η NN – NNTT – TT λ − η η η − η NN – TTTT – T T η η η TT – TT NT – NT Figure 1: Schematic breaking of viscosities in bulk fluid (left) intothe surface viscosities near a flat wall (right). Given are the three/sixnon-zero values found in the components of the stress tensor correla-tor A ijkl , indexed by the normal direction N and the two tangentialdirections T and T (cid:48) . (B) Another symmetry is isotropy , that is the invari-ance of the correlator under rotations. It is not enforcedby the microscopic definition and must result from the av-eraging either over an ensemble or possibly over a longtrajectory. In bulk, the rotations are arbitrary, while ina geometry with plates one direction is distinct from thetwo others and only rotations parallel to the plates are per-mitted. The maximal number of different components isnow reduced to 21. They are, however, not all completelyindependent, and some of them vanish: A systematic de-composition of the rank-four tensor into parallel and or-thogonal parts, using parallel isotropy and Eqs. (7), revealsthat these 21 entries can take only nine different non-zerovalues, and that they depend on seven independent pa-rameters (viscosities).(C) The next symmetry is the invariance under indexpair exchange , (cid:104) s ij s kl (cid:105) = (cid:104) s kl s ij (cid:105) . (8)It is a consequence of the entropy production being ascalar [8]. This symmetry leads to the identification oftwo of the viscosities and the vanishing of another. Afterhaving used all symmetry requirements, we thus have atotal of five viscosities left, which give rise to six differentnon-zero values in the 13 principally different componentsof the stress correlator.The symmetry conditions near a plate require remark-ably more viscosities than those in bulk, where two vis-cosities lead to three different values. The presence of thewall thus induces a separation of the bulk viscosities intosurface viscosities. This splitting is depicted in Fig. 1. Thediagram provides on the left-hand side the three (in bulk)different values found in the components of Eq. (6), andon the right-hand side the six non-zero values near a flatrigid immobile plane. The full rank-four tensor expressingthe random stress correlations near a wall is A ijkl = η N i N j N k N l + η g ij g kl + η ( g ik g jl + g il g jk − g ij g kl ) + η ( N i N j g kl + N k N l g ij )+ η ( N i N k g jl + N j N k g il + N i N l g k + N j N l g ik ) , (9)with the shortcut g ij := δ ij − N i N j . The rightmost column of Fig. 1 uses a symbolic nota-tion for the indices of tensor components. N stands forthe direction normal to the wall, T may stand for bothtangent directions, and T (cid:48) is only used if both tangentdirections are used: In particular, TT stands for the twoindex combinations with the same tangent vector, whereas TT (cid:48) stands for the ones with different tangent vectors.The same notation will be used in the figures of Sec. 4to visualize the 13 components of the stress correlationtensor, six of which are non-zero and seven vanish.
3. Coarse-graining of a microscopic stress tensor
Let us consider a molecular dynamics simulation in or-der to study the microscopic dynamics of the molecules ina liquid. The liquid here consists of point-like atoms, thedynamics of which is governed by Hamilton’s equationscomprising a pair-interaction potential. In this dynamicswe define a microscopic stress tensor . It is chosen suchthat its coarse-grained version , which is obtained by in-tegrating it over a small region, will be consistent with amacroscopic notion of the hydrodynamic stress tensor.The Hamiltonian of N point-like atoms ( α ) of masses m ( α ) , positions r ( α ) ( t ) and linear momenta p ( α ) ( t ), whichundergo mutual influence via a pair-potential φ , reads H = N (cid:88) α =1 p ( α ) · p ( α ) m ( α ) + N (cid:88) α =1 α − (cid:88) β =1 φ (cid:0) (cid:107) r ( β ) − r ( α ) (cid:107) (cid:1) . (10)The time evolutions of the positions and momenta aregiven by Hamilton’s equations. The microscopic stress ten-sor τ ik , can then be introduced as the current density ofthe density of linear momentum g , reading g i ( r , t ) := (cid:88) α p ( α ) i ( t ) δ (cid:0) r − r ( α ) ( t ) (cid:1) . (11)In the absence of external forces, the conservation of mo-mentum is expressed by the following balance equation,which introduces the microscopic stress tensor τ ik , ∂ t g i ( r , t ) − ∂ r k τ ik ( r , t ) = 0 . (12)The definition of these tensor fields allows us to pass froma description in terms of the atom positions r ( α ) ( t ) to afield description for any position r . Equation (12), which isat the same time a conservation theorem and a definition,fixes τ only partially. We may always add an arbitrary partwith zero divergence, a gauge freedom which is inherent tothe continuum description. However, we will not go intothe implications of this freedom. A convenient definitionof the microscopic stress tensor is the following, since itis symmetric and lends itself to a straight-forward coarse-3raining technique, τ ik ( r , t ) = − (cid:88) α p ( α ) i p ( α ) k m ( α ) δ (cid:0) r − r ( α ) ( t ) (cid:1) + 12 (cid:88) α (cid:88) β (cid:54) = α φ (cid:48) (cid:0) (cid:107) r ( β ) − r ( α ) (cid:107) (cid:1) r ( α ) i − r ( β ) i (cid:107) r ( β ) − r ( α ) (cid:107) × (cid:0) r ( α ) k − r ( β ) k (cid:1) (cid:90) δ (cid:16) r − r ( α ) − λ (cid:0) r ( β ) − r ( α ) (cid:1)(cid:17) dλ. (13)The first term is usually called the kinetic part of the stresstensor. The second term comprises virial contributions, to-gether with the integral of a three-dimensional delta func-tion over the line starting at the point r ( α ) and endingat r ( β ) . Expression (13) is the same as that of Forster [16,see Eq. (4.6)], with the exception of the lower boundary ofthis integral.We choose to do a coarse-graining by averaging themicroscopic stress tensor within a small cube V a ( x ) ofside-length a , centred around x . We call this average the coarse-grained stress tensor τ ( a ) ik ( x , t ) at the point x , τ ( a ) ik ( x , t ) := 1 a (cid:90) V a ( x ) τ ik ( r , t ) d r . (14)The integration of the kinetic term in τ ik is readily done. Ituses the momenta of the particles which are found withinthe volume at time t . The line integral of the delta functionhas a simple geometrical interpretation. When integratedover the volume V a ( x ), only those parts of the line integralcontribute where the line is found within the volume. Forthe (convex) cubic box V a ( x ), only four different cases arepossible: (cid:107) r ( β ) − r ( α ) (cid:107) (cid:90) V (cid:90) δ (cid:16) r − r ( α ) − λ ( r ( β ) − r ( α ) ) (cid:17) dλ d r = (cid:107) r ( β ) − r ( α ) (cid:107) if r ( α ) ∈ V and r ( β ) ∈ V , (cid:107) ¯ r − r ( α ) (cid:107) if r ( α ) ∈ V and r ( β ) / ∈ V , (cid:107) r ( β ) − ¯ r (cid:107) if r ( α ) / ∈ V and r ( β ) ∈ V , (cid:107) ¯ r − ¯ r (cid:48) (cid:107) if r ( α ) / ∈ V and r ( β ) / ∈ V . (15)Here, ¯ r and ¯ r (cid:48) denote the points where the connectionline between the atoms intersect the boundary of the vol-ume V a ( x ). The last case in Eq. (15) yields no contri-bution if these intersection points do not exist. Figure 2provides a graphical representation of these different cases.For calculating the coarse-grained stress tensor we thushave to find the contributing atom pairs, to evaluate thepair potential, and to find the intersection points. Noticethat even a box which does not contain any particle mayyield a non-zero microscopic stress. It is sufficient that anyconnecting line passes through the box.The usual way to calculate a stress tensor in molecu-lar dynamics simulations is to add up the kinetic contri-bution and the full virial of the atoms. The simulation ¯ r ¯ r r ( β ) r ( α ) ¯ r V a ( x ) a Figure 2: Visualisation of the virial contributions to the coarse-graining of the stress tensor. From the connection lines betweenthe atoms (small circles) only the bold parts are taken into account.All cases of Eq. (15) are displayed. package
LAMMPS [17], for example, outputs the followingquantity as stress per atom , T ( α ) ik := − m ( α ) v ( α ) i v ( α ) k + 12 (cid:88) β ∈ Nbrs.( α ) ( r ( α ) k − r ( β ) k ) × r ( α ) i − r ( β ) i (cid:107) r ( α ) − r ( β ) (cid:107) φ (cid:48) (cid:0) (cid:107) r ( α ) − r ( β ) (cid:107) (cid:1) . (16)However, T ( α ) ik has the dimension of stress times volume.The stress is attributed to the atom with index α , its neigh-bours have indices β . It is clear that T ( α ) ik cannot be alocal quantity, as it does not depend on any space vari-able. It may be the integral of a local quantity—such asthe stress in Eq. (13)—over a specific volume around theatom α —but over which volume? We find the answer tothis question from the explicit averaging procedure above:Only if the particle and all its neighbours are within thevolume of integration, that is the first case in Eq. (15),then we obtain expression (16) as the contribution com-ing from this atom. We conclude that T ( α ) ik is the bulkcontribution to a stress tensor when integrated over a vol-ume which is sufficiently large that it contains the atomin question and all its interaction partners. It discards allboundary terms. For geometrical reasons it is not possi-ble to assign space-filling individual volumes to all atoms(such as is done for example by a Voronoi tessellation) inorder to obtain the stress per atom as an integral over thelocal stress. As the cell boundaries in the Voronoi tessella-tion are determined only by nearest neighbour’s positions,they cannot correctly represent all interaction partner’spositions, which are spanned over a larger region.As our interest lies in the local fluctuations of the stresstensor, we have to keep the integration boxes as small aspossible. The small boxes forbid the use of Eq. (16), wehave to average the microscopic stress tensor explicitly,according to the cases in Eq. (15). The additional work tobe done in comparison with Eq. (16) is the solution of thegeometrical intersection problem. Throughout the paper, terms specific to
LAMMPS are typesetin a typewriter font. NT (a) T T (b) Figure 3: (Colour on-line) Snapshots of the simulation domain with N = 1380 fluid atoms (black), confined between two fixed plateswhich consist of a regular grid of immobile atoms (red/blue). Shownare cuts through the domain as seen from two different sides. Theseries of cubes in the middle of panel (a) indicate the cubic boxesused for averaging in Sec. 4.2. In the hatched region, a Langevinthermostat controls the temperature.
4. Results from molecular dynamics simulations
The molecular dynamics simulations described in thefollowing have been done with the molecular dynamicspackage
LAMMPS [17]. The coarse-graining procedure de-scribed above has been added as a module.
Our simulation is set up as a three-dimensional Lennard–Jones liquid between two plates, as depicted in Fig. 3.1380 fluid atoms are confined between two fixed plates.The plates are constructed as a single regular layer ofatoms with the same interaction potential. In order toapproximate a smooth flat surface, the lateral distancebetween these plate atoms is about a third of that in thefluid. The simulation domain is periodically continued inthe two directions parallel to the plates.The fluid atoms are advanced in time using the Verletalgorithm (
NVE integration ) which is consistent with themicrocanonical ensemble. An exception are the hatchedregions in Fig. 3a, where an additional Langevin thermo-stat is applied. Such a thermostat is necessary to counter-act numerical integration errors, but it modifies the statis-tics of the stress. We therefore chose to apply it only inregions where we do not measure the stress. The plate atoms are not advanced in time, simply by not applyingany integration routine. The forces on these atoms arecalculated, however, which is equivalent to the possible in-clusion of strong constraining forces from crystalline layersbehind. This procedure allows to determine the force ex-erted on the plates by the fluid.The pair-interaction potential between the atoms is the12 / φ ( r ) = (cid:15) (cid:2) ( σ/r ) − ( σ/r ) (cid:3) for r < . σ r > . σ (cid:80) n = − ξ n ( σ/r ) n else. (17)The six parameters ξ n are determined to guarantee con-tinuous derivatives everywhere, up to the second. We usethis smooth potential in order to make sure that discon-tinuities in the second derivative do not produce artefactsin the measurement of the stress tensor.In the following, we will use Lennard-Jones units forall quantities. In these units, distances are measured asmultiples of the length parameter σ in the Lennard–Jonespotential, energy as multiples of (cid:15) , and time in terms of (cid:112) mσ /(cid:15) , where m is the mass unit. In these units, themass of the atoms is chosen to be 1, the side-length ofthe simulation domain is fixed to be 12. The time unitis approximately the oscillation period in the minimum ofthe pair potential φ . The time-step of integration, used inthe Verlet algorithm is 0 . T = 1 . ρ = 0 . p = 1 .
3. Near the plates,the enforced flat geometry imposes a partial ordering ofthe fluid. Up to three distinct layers of fluid atoms canbe identified near each wall. The first layer has a locallyhexagonal structure as seen in Fig. 3b. Its width and ori-entation are independent of the grid spacing of the plateatoms. Due to periodic boundary conditions, the hexag-onal grid does not rotate anymore as soon as the systemis equilibrated. Of course, also the stress tensor becomesanisotropic close to the plates.The simulation was first run to let the system find itsthermodynamic equilibrium. We then continued the simu-lation while calculating the coarse-grained stress tensors atsome fixed positions every 5 time steps. These time seriesare analysed in the following, either using the blockingmethod [21] to ensure that all used variables are equili-brated, or by simple averaging to measure spatial correla-tions of the stress tensor, or by fast Fourier transforms for5 . r D ˜ τ ( a ) i j ( x , t ) E i n L J - un it s . . . a in LJ-units NN (1 × ) TT (2 × ) TT (1 × ) NT (2 × ) a − / Figure 4: (Colour on-line) Fluctuations of the stress tensor in bulk, asa function of the side a of the averaging box. All tensor componentsare plotted. They are grouped according to the expected symmetry(see text).0 . . . . . . . e s ti m a t e d s t a nd a r dd e v i a ti on o f t h ea v e r a g ee s ti m a t e number of data points in the “blocked” list NN (1 × ) TT (2 × ) TT (1 × ) NT (2 × ) Figure 5: (Colour on-line) Estimated error of the average estimationfrom the blocking transformation [21]. We used the time series ofall stress components in a box of size a = 1 . the time correlations. The averages (cid:104)·(cid:105) in the following arethus time-averages over a long trajectory. As the systemis in equilibrium, this should be equivalent to an ensembleaverage. Before presenting the temporal and spatial dependenceof the stress correlator, let us investigate how the sidelength a of the averaging box changes the fluctuations ofthe coarse-grained stress tensor inside.The explicit coarse-graining method by averaging oversmall cubes allows to observe the transition from thermo-dynamic behaviour (in large boxes) to highly fluctuatingrandom behaviour (in small boxes). The amount of fluc-tuations depends of course on the size of the box. As wetreat here only thermodynamically equilibrated fluids, themean square deviation of the stress tensor should scale asthe inverse volume of the averaging box – just as for anyintensive quantity [22, § (cid:113)(cid:10) ˜ τ ( a ) ij ( x , t ) (cid:11) ∝ a − / , (18)where ˜ τ ( a ) ij ( x , t ) is defined by subtracting the macroscopicstress tensor, which is − p δ ij ,˜ τ ( a ) ij ( x , t ) := τ ( a ) ij ( x , t ) + p δ ij . (19)The value for the pressure is taken from disordered fluidatoms only, and not from the layers close to the plates.The numerical verification of the scaling (18) is plottedin Fig. 4. Amazingly, we find the correct scaling already forbox sizes of a ∼
2, that is a side length of the order of theaverage distance between the fluid atoms. The scaling (18)is originally based on a thermodynamic argument whichrequires that many particles are involved.
Instead, we heresee the correct scaling already for very few particles . Itappears that there is more truth in the thermodynamicargument than just the interaction of many particles.For visualising the tensor of stress fluctuations, we havechosen the following strategy, which will be applied also tothe following plots: The symmetry requirements (B) and(C) in Sec. 2 lead to components which are equal only in astatistical sense. All components which may exhibit suchfluctuations are plotted by individual lines. The legendand the colour coding, however, differentiate only betweenthe groups of components which should be equal on av-erage. The number of lines actually plotted is given inparentheses.The fact that there is no difference between the nor-mal and the parallel directions in the data of Fig. 4 con-firms our assumption that the fluid behaviour in the mid-dle between the walls is the same as if they were absent.There, the diagonal and the off-diagonal parts fluctuatedifferently. This is a consequence of the components beingdifferent functions of the two viscosities (see also Fig. 1).Exemplary for all plotted variables, we show for thecentred boxes that the system is well equilibrated. Fig. 5visualises the blocking transformation [21] of the diagonalentries in the coarse-grained stress tensor. A plateau isclearly visible, which shows that the analysed time seriesindeed correspond to a stationary stochastic process.
The aim of this paper is to determine the spatial andthe temporal structure in correlations of fluctuations ofthe coarse-grained stress tensor ˜ τ ( a ) ik . We are interested inthe difference between correlations in bulk and in the nearvicinity of the walls. In order to do so, we set out threedifferent series of cubic boxes, as depicted in Fig. 6. Thesecubes are shifted by small distances and may overlap. Thefirst series is in the middle between the plates, where thefluid behaves as in bulk, a second one close to a wall andparallel to it, and a third series is normal to the wall.The displacement of the boxes makes it possible todetermine the spatial form of the stress correlation by6 igure 6: (Colour on-line) Same as Fig. 3a, but with different seriesof overlapping coarse-graining cubes, which are used in Figs. 7–10. cross-correlating the coarse-grained stress fluctuations indifferent boxes. Temporal correlation functions will alsobe given. We would like to compare the result with theuncorrelated theory in Eq. (5). However, the nonzero sizeof the boxes used for averaging does not allow a directcomparison. Remind that ˜ τ ( a ) ij depends on the box size.We should find the delta function of Eq. (5) in the limit a →
0, which is infeasible in practice. For all nonzerosizes, the numerical data is sort of “convoluted” with theshape of the boxes. The best we can do is to compare withthe equally smoothed version of Eq. (5) using a reasonablysmall a . The corresponding correlator of the uncorrelatedstress tensor s ij is expressed by the relative overlap of thetwo boxes, (cid:10) ˜ τ ( a ) ij ˜ τ ( a ) kl (cid:11) ∝ (cid:12)(cid:12) V a ( x ) ∪ V a ( y ) (cid:12)(cid:12) , (20)of course disregarding the infinite term δ ( t − s ) in Eq. (5).Here, | · | denotes taking the volume of a set.The spatial and the temporal correlations of the ran-dom coarse-grained stress tensor are depicted in Fig. 7.The data have been generated by correlating the contentseach box in the central series of Fig. 6 with the first one,which has x as its centre. Plotted are all 36 componentsof the resulting rank-four tensor which can be different us-ing the symmetry (A) of Sec. 2. They are grouped intothe 13 different classes according to the symmetries (B)and (C). Fig. 7 shows that seven of them vanish (greylines) and that the remaining six take three different val-ues. This is expected according to what has been saidaround Fig. 1.The spatial correlations in Fig. 7a should be comparedwith the results expected from an uncorrelated stress ten-sor. These would be lines starting at some value whichis related to the viscosities (unspecified due to δ ( t − s ) inEq. (5)) and which goes to zero the distance a and beyond.The coincidence in Fig. 7a is remarkable. The numericvalues are a bit more smoothed, but the apparent convo-lution kernel responsible for this smoothing has a width ofless than unity. We therefore conclude that Eqs. (5) withEq. (6) in bulk are very precise, down to a length scale of the order of the molecule distance.
Similarly to what wefound in Sec. 4.2, such a strong result cannot be expectedfrom thermodynamic considerations.The temporal correlations in Fig. 7b exhibit the devi-ations from the temporal delta function in Eq. (5). Thedecay is sufficiently fast to reasonably adopt the approx-imation of the delta function. Whether there is an alge-braic contribution, the so-called long-time tail cannot bejudged from these data. In order to calculate the viscosityvalues, we use the integral over the detailed time correla-tor, as proposed by the Green–Kubo equation. The ideais the same as in Refs. [23, 24, 25], only that we take a boxsmaller than the simulation box. For the shear viscosity,this reads η = a kT ∞ (cid:90) dt (cid:68) ˜ τ ( a ) NT ( x ,
0) ˜ τ ( a ) NT ( x , t ) (cid:69) , (21)and similarly for the two components comprising the vol-ume viscosity λ . From the data in Fig. 7b we find thevalues η = 0 .
84 and λ = 1 .
72 by regression. Near the plates, the stress correlations look quantita-tively and qualitatively different. The average stress ten-sor and also its fluctuations are anisotropic. Our aim inthis section is to determine the five surface viscosities inorder to quantify their deviations from the bulk values,according to the scheme in Fig. 1.First, let us look at the coarse-grained random stresstensor itself. Fig. 8 reveals the mechanical properties ofthe first fluid layers: We find a huge lateral diagonal stress( TT components) which strongly depends on the distanceto the wall but not on the lateral position. In the fig-ure, the points x correspond to the centres of the seriesof boxes laid out in Fig. 6. The position x is the centreof the first box in the normal series, which is closest tothe wall. The average value of the huge diagonal stress iseven larger than its fluctuations. This lateral stress canbe understood as a surface tension of the contact betweenliquid and solid. The NN -component is much smaller—itsaverage value is the same as the bulk pressure. The NN -component is dominated by its fluctuations which grow inthe vicinity of the wall and even overgrow the tangentialstresses. In the following, we will find that there are cor-relations in these fluctuations. The lines in Fig. 8b provethat both the average values and the standard deviationsare independent of the lateral position. The roughness ofthe NN -component in this plot is due to the large fluctu-ations. Notice that these values do not quite correspond to the findingin Refs. [24, 25]. A possible explanation is the dependence of theviscosities on the size a of the averaging box. We found such adependence together with long-time tails in the autocorrelation ofthe random stress tensor in bulk. D ˜ τ ( a ) i j ( x , t ) ˜ τ ( a ) k l ( x , t ) E i n L J - un it s . . . k x − x k in LJ-units NN – NN (1 × ) TT – TT (2 × ) TT – T T (2 × ) NN – TT (4 × ) TT – TT (1 × ) NT – NT (2 × ) NN – NT (4 × ) TT – TN (4 × ) TT – TT (4 × ) NN – TT (2 × ) TT – T N (4 × ) NT – NT (2 × ) TN – TT (4 × ) Σ (cid:0) ˜ τ ( a ) i j ˜ τ ( a ) k l (cid:1) (a) − . . . . . D ˜ τ ( a ) i j ( x , t ) ˜ τ ( a ) k l ( x , s ) E i n L J - un it s . . . . . . . . . | t − s | in LJ-units NN – NN (1 × ) TT – TT (2 × ) TT – T T (2 × ) NN – TT (4 × ) TT – TT (1 × ) NT – NT (2 × ) NN – NT (4 × ) TT – TN (4 × ) TT – TT (4 × ) NN – TT (2 × ) TT – T N (4 × ) NT – NT (2 × ) TN – TT (4 × ) Σ (cid:0) ˜ τ ( a ) i j ˜ τ ( a ) k l (cid:1) .
25 0 . .
75 1 (b)
Figure 7: (Colour on-line) Spatial and temporal correlations of the stress fluctuations in the middle of the domain (in bulk). Plotted are theaverage correlators and their standard deviations (Σ( · ), see insets) for different classes of index combinations. For the data in Fig. 8 and in the following ones, anothercomment on the symmetry is necessary. By the noted sym-metries in Sec. 2, the two lines for the TT -components inFig. 8 should coincide. However, they do not. The reasonis that the structured liquid is not entirely isotropic in thetangential directions, but it is ordered in a hexagonal grid,as depicted in Fig. 3b. This strict ordering is a result ofthe wall being very flat. A roughness of the size of oneatom distance would easily destroy the hexagonal struc-ture, and only the symmetries of Sec. 2 would remain. Inorder to weaken the effect, we averaged over eight trajecto-ries which differed only by the orientation of the hexagonalgrid. This reduced the deviations considerably to those inFig. 8.The spatial correlations of the random stress tensorare depicted in Fig. 9. Again, the two series of integra-tion boxes in different directions are used, see Fig. 6. Theplot shows the correlation of each cube in the series withthe first one, which is in the normal series the one closelyaligned at the wall. For the normal series, the result is notexpected to be independent of this choice. Only in the par-allel series, the correlation should depend only on the dis-tance of the boxes. As explained above, there are in prin-ciple 13 different results to be expected. Seven of them,those indicated in grey, are zero. The six non-zero compo-nents are plotted in colour. The correlations of the stresstensor inherit the main features of the stress in Fig. 8, inparticular the systematic part which we called a surface tension. In normal direction, the tangential componentsexhibit a similar anti-correlation as in Fig. 8b. The decayof the most pronounced correlations takes place on thelength of one atom-distance and is independent of the boxsize a . We verified this statement using other box sizes. Intangential direction, the surface tension manifests itself ina nonzero value of the tangential correlation, see Fig. 8b.Apart from the surface tension, we find in Fig. 8b asimilar behaviour as in bulk, that is a linear decay up toa distance a and zero beyond. This is clearly visible forthe NN – NN component, but also for the others—aftersubtracting the contribution of the surface tension. Thisimplies that the spatial delta function still works parallelto the walls. The correlation in normal direction, which isdominated by the details of the interaction between walland fluid, depends both on the distance between the boxand their distance to the wall. These two considerationslead us to the following proposition for the correlationfunction of the random stress tensor near an immobile flatrigid wall: (cid:10) s ij ( x , t ) s kl ( y , s ) (cid:11) = 2 kT A ijkl δ ( t − s ) × δ ( x || − y || ) f ( x ⊥ , y ⊥ ) . (22)with the tensor A ijkl from Eq. (9), depending on the sur-face viscosities, and with x || and x ⊥ the parallel and nor-mal projections, respectively. The function f , which hasthe physical dimension of inverse length, takes into ac-8 − D ˜ τ ( a ) i j ( x , t ) E i n L J - un it s . . . k x − x k in LJ-units NN (1 × ) TT (2 × ) TT (1 × ) NT (2 × ) Σ (cid:0) ˜ τ ( a ) i j ( x , t ) (cid:1) (a) − − D ˜ τ ( a ) i j ( x , t ) E i n L J - un it s . . . k x − x k in LJ-units NN (1 × ) TT (2 × ) TT (1 × ) NT (2 × ) Σ (cid:0) ˜ τ ( a ) i j ( x , t ) (cid:1) (b) Figure 8: (Colour on-line) The coarse-grained stress tensor as a func-tion of the position, as indicated in two series close to the wall inFig. 6: (a) parallel to the plate, (b) orthogonal to the plate. count the microscopic details of the interaction betweenwall and fluid. As the only length scale involved is thedistance between wall and the first structured liquid layer,we expected f to scale as one over this distance.It remains to check whether the delta function in timeis still a good approximation, and what values the surfaceviscosities have. The temporal correlations near the wallare depicted in Fig. 10. The results are taken from one ofthe boxes in the parallel series. Here, we see even betterthan in Fig. 9 that there is a systematic contribution com-ing from the surface tension. Notice that it affects onlythe tangential components and not the heavily oscillating NN – NN component, which is decaying to zero, but on alonger time scale than plotted.In order to obtain the viscosities from these data, wehave subtracted the asymptotic value from the correla-tion function, which renders the function integrable. Thismeans that we implicitly subtracted a sort of local pressurewhich is allowed to be anisotropic. This pressure reflectsthe static contributions from the interaction liquid–wall,which we called surface tension above. The Green–Kubointegral of the such normalized correlation functions yields the following values for the surface viscosities: η ≈ ,η between 16 and 34 ,η ≈ . ,η between 65 and 178 ,η between 4.2 and 5.5 . (23)The values for these viscosities are not very precise forthree reasons: The first is the hexagonal structure of thefirst fluid layer which does not coincide with the Cartesiancoordinates used for the tensor indices. We thus expectdifferent viscosity values for different orientations of thehexagonal grid. The above mentioned average over eightdifferent trajectories with different orientations apparentlydoes not suffice to obtain unified values. The second rea-son is that the thermodynamic pressure is not isotropicanymore and that the time-correlation function does notdecay to zero. We thus had to strip off a linear time-dependency before integrating, as already explained. Thethird reason is that the intrinsic time scales of the time-correlation function are much longer than in bulk. Whilethe bulk viscosities could be easily identified from a Green–Kubo integral over a few dozens of timesteps, it requiredthousand and more timesteps to do the same for the sur-face viscosities. This much longer time scale can alreadybe guessed from the comparison of Figs. 7b and 10. The rigid plates used in the above simulations consist ofthe same type of atoms as the liquid. Such walls are knownto be hydrophilic, which is also seen in the strict ordering ofthe first structured fluid layer. We repeated the simulationusing hydrophobic walls, which have a modified Lennard–Jones potential to interact with the fluid atoms [26], φ ( r ) = 4 (cid:15) (cid:2) ( σ/r ) − c ( σ/r ) (cid:3) (24)with the same cutoff and with an equivalent smoothingof the potential as in Eq. (17). According to Barrat andBocquet [26], who proposed the form (24) for simulatinghydrophobic walls, the factor factor c = 0 . ◦ .The simulation leads to a similarly structured liquid asin Fig. (3), only that the first layers are a bit farther awayfrom the walls, that they are not as flat, and that their lat-eral structure is not as regular. The ordering effect of thefirst layers are nevertheless very strong, and must there-fore be induced by the flat geometry of the wall. We tookthe same bulk pressure and density as in the previous sim-ulation. As expected for a hydrophobic surface, the TT -components in Fig. 8 are different. They are considerablysmaller, around the value 3 instead of 13. The correlationfunctions look similar to those depicted in Figs. 9 and 10,with the difference that they decay much faster in time,in particular the NN – NN component. Accordingly, also9 D ˜ τ ( a ) i j ( x , t ) ˜ τ ( a ) k l ( x , t ) E i n L J - un it s . . . k x − x k in LJ-units NN – NN (1 × ) TT – TT (2 × ) TT – T T (2 × ) NN – TT (4 × ) TT – TT (1 × ) NT – NT (2 × ) NN – NT (4 × ) TT – TN (4 × ) TT – TT (4 × ) NN – TT (2 × ) TT – T N (4 × ) NT – NT (2 × ) TN – TT (4 × ) Σ (cid:0) ˜ τ ( a ) i j ˜ τ ( a ) k l (cid:1) (a) − D ˜ τ ( a ) i j ( x , t ) ˜ τ ( a ) k l ( x , t ) E i n L J - un it s . . . k x − x k in LJ-units NN – NN (1 × ) TT – TT (2 × ) TT – T T (2 × ) NN – TT (4 × ) TT – TT (1 × ) NT – NT (2 × ) NN – NT (4 × ) TT – TN (4 × ) TT – TT (4 × ) NN – TT (2 × ) TT – T N (4 × ) NT – NT (2 × ) TN – TT (4 × ) Σ (cid:0) ˜ τ ( a ) i j ˜ τ ( a ) k l (cid:1) (b) Figure 9: (Colour on-line) Spatial correlations of the random stress tensor in the vicinity of a wall. (a) in normal direction, and (b) parallelto the wall, see the series of cubes in Fig. 6. Each box in the series is correlated with the first one, centred around x (in the orthogonal case:the one closest to the wall). Legend and inset as in Fig. 7a. the surface viscosities turn out to be different. Near thehydrophobic surfaces, we find the values η ≈ ,η ≈ . ,η ≈ . ,η ≈ . ,η ≈ . . (25)It is not astonishing that these viscosities are smaller thanthose found near the wetting walls. In particular in tan-gential direction, where the liquid is feels less constraints,the dissipation is smaller. This affects η – η , but not η , which is accordingly of a similar magnitude as in thewetting-wall case. Particularly astonishing is the smallvalue of η , which takes over the role of the shear viscos-ity. It is much smaller than the bulk value ( η = 0 .
84 alsoin the hydrophobic simulation).
5. Discussion
The numerical observations above, as described in Sec-tions 4.3 and 4.4 allow to estimate the validity of the mainassumptions leading to fluctuating hydrodynamics. In-deed, we observed that the spatial delta function in Eq. (5) works very well in bulk. Near a plate, which has beenassumed to be flat, immobile and rigid, the isotropy isrestricted to rotations parallel to the wall, and we mustreplace the spatial delta function by one which is only act-ing parallel to the wall, as given in Eq. (22). In orthogonaldirection, the details of the dynamical interaction betweenliquid and wall introduce an unknown factor in form ofa function depending on the orthogonal distances. More-over, we have to deal with five viscosities instead of two,the values of which are provided in Eq. (23).The parallel structure of the spatial delta functions andthe number of viscosities are both imposed by pure sym-metry and cannot be disputed. But we may ask whichof the surface viscosities will have an effect on measur-able quantities, such as on the flow profile in nanoporesor on the motion of a Brownian particle. If we regardonly forces exerted on the rigid walls—and their correla-tion functions—, there are only two viscosities which mayplay a role, namely those possessing a normal direction inboth contributing stresses. They are easily identified inthe scheme in Fig. 1 to be η and η , the latter of whichtakes over the role of the bulk shear viscosity, and thefirst one the role of the bulk volume viscosity. It is in-teresting to see that the value of η is entirely indepen-dent on the shear viscosity in bulk. It depends more onthe hydrophobic/hydrophilic character of the walls. The10 D ˜ τ ( a ) i j ( x , t ) ˜ τ ( a ) k l ( x , s ) E i n L J - un it s . . . . . . . . . | t − s | in LJ-units NN – NN (1 × ) TT – TT (2 × ) TT – T T (2 × ) NN – TT (4 × ) TT – TT (1 × ) NT – NT (2 × ) NN – NT (4 × ) TT – TN (4 × ) TT – TT (4 × ) NN – TT (2 × ) TT – T N (4 × ) NT – NT (2 × ) TN – TT (4 × ) Σ (cid:0) ˜ τ ( a ) i j ˜ τ ( a ) k l (cid:1) .
25 0 . .
75 1Figure 10: (Colour on-line) Temporal correlations of the random stress tensor near a wall. Legend and inset as in Fig. 7b. surface-equivalent of the volume viscosity, η , is in bothcases nearly two orders of magnitude larger than in bulk.The finding that the surface viscosities do not havevalues comparable to the bulk viscosities brings us nowto more profound questions on the dynamical theory in-volved: For the Brownian particle, we argue that only thefirst surrounding liquid layers do actually interact with theparticle. There, the relevant dissipation coefficients are notthe bulk viscosities but the two viscosities η and η . Itappears thus as a reasonable to ask why we do not see η instead of η appearing in the Stokes–Einstein formula forthe diffusion coefficient? We have to leave this questionopen at the moment. A possible answer might go in thesame direction as the discussion of time correlations andspatial ordering effects of Ref. [27].A possible answer is certainly linked to the fact thatthe surface viscosities here are local quantities. They havebeen determined using coarse-graining over very small boxes.The usual definition of a “viscosity”, however, implies lim-its to large boxes and long time integrals (or the limit ω → k → in Fourier space and time [3,p. 186]). In this limit, the five surface viscosities cannot beidentified anymore. Another implication of the use of lo-cal viscosities , which are not viscosities in the usual senseof a large-box limit, is that they do not lead to differ-ential equations. They rather appear as integral kernels,nonlocal both in space and in time, in integro-differentialequation. We cannot say much about the possibility tofind an to interpret such equations at the moment, theless as we expect the long-time tails in the autocorrelationfunction, thus nonlinear terms, to intervene. Neverthe-less, we would like to stress that we do not expect thequalitative nature of the anisotropy described here—thatis the large difference of surface and bulk viscosity values,and that the surface viscosities do not depend on the bulkviscosities—to depend on the precise temporal behaviour.The long-time tails appearing in the correlation functionsshould only change the quantitative values of the viscosi-ties, leaving the presented material untouched. Acknowledgements
I would like to dedicate this paper to Peter H¨anggi,on occasion of his 60th birthday, which is celebrated bythe present special issue of Chemical Physics. The subjectof this paper, in particular the aspects explained in thelast section touches in part his early works on stochasticprocesses with memory.I would further like to express my thanks to PeterTalkner for waking my interest in fluctuating hydrodynam-ics, and to Laurent Joly and Lyd´eric Bocquet for intro-ducing me to the molecular dynamics package
LAMMPS .Further sincere thanks go to Anthony C. Maggs for point-ing me to the blocking method and to Claire Lemarchand,Ken Sekimoto, Falko Ziebert and J¨org Rottler for helpfuldiscussions.
References [1] L. D. Landau, E. M. Lifshitz, Fluid Mechanics, Pergamon Press,1959.[2] R. F. Fox, G. E. Uhlenbeck, Contributions to Non-EquilibriumThermodynamics. I. Theory of Hydrodynamical Fluctuations,1970 13 (Phys. Fluids.) 1893.[3] R. Kubo, M. Toda, N. Hashitume, Statistical Physics II.Nonequilibrium Statistical Mechanics, Springer Verlag, 1985.[4] H. Grabert, P. Talkner, P. H¨anggi, Microdynamics and Time-Evolution of Macroscopic Non-Markovian Systems, Z. Physik B26 (1977) 389.[5] E. H. Hauge, A. Martin-L¨of, Fluctuating hydrodynamics andBrownian motion, J. Stat. Phys. 7 (1973) 259.[6] T. S. Chow, J. J. Hermans, Effect of Inertia on the BrownianMotion of Rigid Particles in a Viscous Fluid, J. Chem. Phys. 56(1972) 3150.[7] G. Temple, Cartesian Tensors. An Introduction, Methuen, Lon-don, 1960.[8] S. R. de Groot, P. Mazur, Non-Equilibrium Thermodynamics,Dover Publications, New York, 1984.[9] C. Kettner, P. Reimann, P. H¨anggi, F. M¨uller, Drift Ratchet,Phys. Rev. E 61 (2000) 312.[10] Z. Guttenberg, A. Rathgeber, S. Keller, J. O. R¨adler, A. Wix-forth, M. Kostur, M. Schindler, P. Talkner, Flow profiling ofa surface-acoustic-wave nanopump, Phys. Rev. E 70 (2004)056311.
11] M. Schindler, P. Talkner, M. Kostur, P. H¨anggi, Accumulatingparticles at the boundaries of a laminar flow, Physica A 385(2007) 46.[12] K. John, P. H¨anggi, U. Thiele, Ratched-driven fluid transportin bounded two-layer films of immiscible liquids, Soft Matter 4(2008) 1183.[13] L. Bocquet, E. Charlaix, Nanofluidics, from bulk to interfaces,Chem. Soc. Rev. 39 (2010) 1073.[14] M. Schindler, Boundary-integral formulation of particle motionin unsteady Stokes/Oseen flow, submitted .[15] S. Bohlius, H. R. Brand, H. Pleiner, Macroscopic dynamics ofuniaxial magnetic gels, Phys. Rev. E 70 (2004) 061411.[16] D. Forster, Hydrodynamic fluctuations, broken symmetry, andcorrelation functions, Addison-Wesley, 1990.[17] S. J. Plimpton, Fast Parallel Algorithms for Short-Range Molec-ular Dynamics, J. Comp. Phys. 117 (1995) 1.[18] B. Smit, Phase diagrams of Lennard–Jones fluids, J. Chem.Phys. 96 (1992) 8639.[19] M. A. van der Hoef, Free energy of the Lennard–Jones solid,J. Chem. Phys. 113 (2000) 8142.[20] E. A. Mastny, J. J. de Pablo, Melting line of the Lennard–Jones system, infinite size, and full potential, J. Chem. Phys.127 (2007) 104504.[21] H. Flyvbjerg, H. G. Petersen, Error estimates on averages ofcorrelated data, J. Chem. Phys. 91 (1989) 461.[22] L. D. Landau, E. M. Lifshitz, Statistical Physics, Part 1,Butterworth-Heinemann, Oxford, 3rd edn., 2001.[23] R. L. Rowley, M. M. Painter, Diffusion and Viscosity Equationsof State for a Lennard–Jones fluid Obtained from MolecularDynamics Simulations, Int. J. Thermophys. 18 (1997) 1109.[24] K. Meier, A. Laesecke, S. Kabelac, Transport coefficients of theLennard–Jones model fluid. I. Viscosity, J. Chem. Phys. 121(2004) 3671.[25] K. Meier, A. Laesecke, S. Kabelac, Transport coefficients of theLennard–Jones model fluid. III. Bulk viscosity, J. Chem. Phys.122 (2005) 014513.[26] J.-L. Barrat, L. Bocquet, Large Slip Effect at a NonwettingFluid–Solid Interface, Phys. Rev. Lett. 82 (1999) 4671.[27] A. V. Mokshin, R. M. Yulmetyev, P. H¨anggi, Diffusion processesand memory effects, New J. Phys. 7 (2005) 9.11] M. Schindler, P. Talkner, M. Kostur, P. H¨anggi, Accumulatingparticles at the boundaries of a laminar flow, Physica A 385(2007) 46.[12] K. John, P. H¨anggi, U. Thiele, Ratched-driven fluid transportin bounded two-layer films of immiscible liquids, Soft Matter 4(2008) 1183.[13] L. Bocquet, E. Charlaix, Nanofluidics, from bulk to interfaces,Chem. Soc. Rev. 39 (2010) 1073.[14] M. Schindler, Boundary-integral formulation of particle motionin unsteady Stokes/Oseen flow, submitted .[15] S. Bohlius, H. R. Brand, H. Pleiner, Macroscopic dynamics ofuniaxial magnetic gels, Phys. Rev. E 70 (2004) 061411.[16] D. Forster, Hydrodynamic fluctuations, broken symmetry, andcorrelation functions, Addison-Wesley, 1990.[17] S. J. Plimpton, Fast Parallel Algorithms for Short-Range Molec-ular Dynamics, J. Comp. Phys. 117 (1995) 1.[18] B. Smit, Phase diagrams of Lennard–Jones fluids, J. Chem.Phys. 96 (1992) 8639.[19] M. A. van der Hoef, Free energy of the Lennard–Jones solid,J. Chem. Phys. 113 (2000) 8142.[20] E. A. Mastny, J. J. de Pablo, Melting line of the Lennard–Jones system, infinite size, and full potential, J. Chem. Phys.127 (2007) 104504.[21] H. Flyvbjerg, H. G. Petersen, Error estimates on averages ofcorrelated data, J. Chem. Phys. 91 (1989) 461.[22] L. D. Landau, E. M. Lifshitz, Statistical Physics, Part 1,Butterworth-Heinemann, Oxford, 3rd edn., 2001.[23] R. L. Rowley, M. M. Painter, Diffusion and Viscosity Equationsof State for a Lennard–Jones fluid Obtained from MolecularDynamics Simulations, Int. J. Thermophys. 18 (1997) 1109.[24] K. Meier, A. Laesecke, S. Kabelac, Transport coefficients of theLennard–Jones model fluid. I. Viscosity, J. Chem. Phys. 121(2004) 3671.[25] K. Meier, A. Laesecke, S. Kabelac, Transport coefficients of theLennard–Jones model fluid. III. Bulk viscosity, J. Chem. Phys.122 (2005) 014513.[26] J.-L. Barrat, L. Bocquet, Large Slip Effect at a NonwettingFluid–Solid Interface, Phys. Rev. Lett. 82 (1999) 4671.[27] A. V. Mokshin, R. M. Yulmetyev, P. H¨anggi, Diffusion processesand memory effects, New J. Phys. 7 (2005) 9.