Shear Thinning in Lennard-Jones Fluids by Stochastic Dissipative Molecular Dynamics Simulation

# Shear Thinning in Lennard-Jones Fluids by Stochastic Dissipative Molecular Dynamics Simulation

Phil Attard
###### Abstract

The shear viscosity of a Lennard-Jones fluid is obtained by stochastic dissipative molecular dynamics (SDMD) simulations. A generic constraint to the equations of motion is given that reduces the sensitivity of the shear viscosity to the value of the the fluctuation-dissipation or thermostat parameter. At high shear rates the shear viscosity becomes dependent on the size of the system, and corrections to the equipartition kinetic temperature arise. At constant kinetic temperature the shear viscosity is shown to decrease with increasing shear rate.

## I Introduction

The value of computer simulations is that they provide molecular-level detail, and quantitative results that are completely reliable. In principle the only approximations invoked are the finite-size of the system and the finite statistical sampling of the configurations, both of which which can be increased to improve or check the accuracy of the results.Allen87 ()

In practice the situation can be more complicated, and it is debatable whether simulations are completely reliable for non-equilibrium systems, since for these cases there is little agreement on the foundations or formulation of statistical mechanics or thermodynamics. Contention continues about fundamental matters such as the form of the probability distribution, the nature of the appropriate equations of motion, and the physical origin and meaning of various thermostats. Whilst in the linear regime different non-equilibrium simulation methods give consistent results, in the non-linear regime there is considerable uncertainty and quantitative disagreement between different computer simulation methods.

For example, the first computer simulations of shear flow in simple atomic fluids were carried out over 30 years ago,Erpenbeck84 (); Heyes86 () but the many papers on the subject since then have produced contradictory results in the high shear rate regime, with both shear thickening, and shear thinning observed (see, for example, Refs. [Todd07, ; Hoover08, ] and references therein), as well as sensitivity to the particular thermostat employed. Hoover08 (); Khare97 (); Berro11 (); Delhommelle03 () Even the structure of shearing simple fluids has been disputed, with some simulations supporting the existence of a string phase at high shear rates, Erpenbeck84 (); Loose89 (); Loose92 (); Bagchi96 () and others showing the opposite. Evans86 (); Padilla95 (); Travis96 (); Delhommelle03 ()

The non-equilibrium molecular dynamics (NEMD) method dominates the field, but in contrast to equilibrium simulations it is arguable (see Ref. [NETDSM, ] and below) that this is based on five explicit approximations of contentious validity: (1) The equations of motion are created artificially to give the desired linear flow, but they have no thermodynamic justification or demonstrated validity in the non-linear regime. (2) The thermostats that have been used likewise have no thermodynamic basis (arguments based on a local Lagrangian co-moving reference frame neglect the effects of any velocity gradient), and they variously invoke a non-physical extra phase space variable, or an artificial constraint on the kinetic energy. (3) Every thermostat relies upon a formula for the local temperature, most usually the equipartition theorem or Rugh’s variant, which formulae are invalid in the non-linear non-equilibrium regime (see Eqs (II.1.2) and (2.33) below). (4) Applying a thermostat assumes that the temperature should be uniform, which need not in fact be the case in a non-equilibrium system (see Eq. (2.36) below). (5) The stochastic contribution to the equations of motion, which necessarily results from the projection of the reservoir interactions, is neglected.NETDSM () The effects of these five approximations are tolerable or negligible in the linear regime, but their unknown consequences in the non-linear regime undermine confidence in the results. Indeed, Hoover, one of the inventors of NEMD, finds that at high shear rates the DOLLS and SLLOD equations of motion for shear flow quantitatively and even qualitatively disagree with each other and with a more physical model of boundary driven flow.Hoover08 ()

This paper presents non-equilibrium computer simulation results for a shearing Lennard-Jones fluid. The concrete aim is to obtain the shear viscosity as a function of shear rate, and to establish as reliably as possible whether or not a simple fluid is shear thinning or shear thickening. More broadly, the analysis and results shed light on some of the reasons for the discrepancies in previous simulations at high shear rates. It is found analytically and numerically that a shearing fluid is inhomogeneous in density and temperature, and consequently that the shear viscosity can be dependent on the system size (as also is the pressure). It is also found that the viscosity can be sensitive to the value of the fluctuation-dissipation thermostat parameter, which sensitivity can be reduced by adding a constraint to the equations of motion. Finally, it is found that whether the fluid is shear thinning or shear thickening depends to some extent on whether the shear flow is boundary driven or else uniformly driven, and on the sign of the temperature derivative of the viscosity.

Of course one can legitimately ask whether the present results are any more reliable in the non-linear regime than previous NEMD and other results. Erpenbeck84 (); Heyes86 (); Todd07 (); Hoover08 (); Khare97 (); Berro11 (); Delhommelle03 (); Loose89 (); Loose92 (); Bagchi96 (); Evans86 (); Padilla95 (); Travis96 () This question is addressed in the conclusion, but here it may simply be pointed out that the present stochastic dissipative molecular dynamics (SDMD) algorithm is a representation of the formally exact transition probability of a non-equilibrium system, and as such it preserves the non-equilibrium probability in phase space. NETDSM (); STD2 (); AttardV () The SDMD method has previously been tested for conductive heat flowAttardV () and for driven Brownian motion.AttardIX () For the present case of shear flow, the SDMD approach is essentially the same as the stochastic Langevin equation, which has previously been used for molecular dynamics simulations of complex shear flow, particularly of bead polymers (see Refs [Kremer90, ; Lyulin99, ; Shang17, ] and references therein). In contrast to the NEMD and similar methods, the SDMD algorithm avoids artificial equations of motion, and unjustified thermostats. Neither does it have to postulate an equilibrium functional form nor a particular local value for the temperature. Finally, and uniquely, it explicitly accounts for the projected reservoir interactions.

In principle the SDMD results are exact, but again the situation in practice is more complicated. Even with the constraint introduced below, large values of the fluctuation-dissipation parameter still affect the results, which creates uncertainty at high shear rates. There is also a certain conceptual or physical uncertainty that arises from the inhomogeneous nature of the shearing system, and from the size-dependence of the viscosity, which lead to questions regarding whether a boundary driven model or a uniform shearing model is most appropriate. These limit confidence in the results at very high shear rates, and are discussed in more detail in the conclusion. Nevertheless, the results at moderate shear rates are sufficiently clear to conclude that departing the linear regime Lennard-Jones fluids are shear thinning.

## Ii Equations of Motion and Simulation Details

### ii.1 General Formulation

In general, for a non-equilibrium system, the stochastic, dissipative equations of motion in phase space are NETDSM (); STD2 (); AttardIX ()

 q(t+Δt) = q(t)+Δtmp(t) p(t+Δt) = p(t)+Δt˙p0(t)+R(t). (2.1)

Here are the positions, with , and are the momenta, with being each particle’s mass. The adiabatic derivative of the momentum is given by Hamilton’s equations, . These equations preserve the non-equilibrium phase space probability density to linear order in the time step , and in the fluctuation-dissipation parameter .

The reservoir force (times the time step) comprises dissipative and stochastic parts, . The dissipative part isSTD2 ()

 ¯¯¯¯¯R(t)=θ|Δt|2kB∇pSr,st(Γ)+θ|Δt|2kB[ˆt−1]¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯∇pSr,st, (2.2)

and the stochastic part is Gaussian distributed with zero mean and variance

 ⟨~R(t′)~R(t)⟩=θ|Δt|δ(t−t′)I––. (2.3)

In these , and is the fluctuation-dissipation parameter. (A constraint may be added to the dissipative force to ameliorate the direct influence of the reservoir on the transport coefficient; see §II.1.1 below.)

The total entropy is the sum of the sub-system entropy and the reservoir entropy. For a point in phase space, the sub-system entropy is zero, .TDSM () For a non-equilibrium system, the reservoir (hence total) entropy comprises static and dynamic parts, . Obviously, the phase space probability for the non-equilibrium system is . NETDSM (); STD2 (); AttardV ()

The static part of the reservoir entropy is the entropy that would be written down if the sub-system were instantaneously an equilibrium system. If the exchangeable linear additive variables are , and the conjugate field variables are , then for the reservoir . In the usual Gibbsian fashion, the static part of the reservoir entropy is

 Sr,st(Γ,t)=−Xs(Γ,t)⋅xr. (2.4)

This procedure is concretely illustrated with shear flow below.

The dynamic part of the the reservoir entropy is NETDSM (); STD2 (); AttardV ()

 Sr,dyn(Γ,t)=−∫t−∞dt′˙S0r,st(¯¯¯¯Γ(t′|Γ,t),t′), (2.5)

where in the integrand appear the adiabatic time derivative of the static reservoir entropy, and also the most likely trajectory from the current point. The dynamic part of the reservoir entropy is not required explicitly in the present paper.

The fluctuation-dissipation parameter in Eqs (2.2) and (2.3) can be regarded as a thermodynamic drag coefficient that occurs in the transition probability. Its value cannot be too large since the stochastic dissipative equations of motion are based upon an expansion to linear order in , but it is otherwise arbitrary. The essential point is that it must link the dissipative and stochastic forces exactly as in Eqs (2.2) and (2.3). This is the fluctuation-dissipation theorem, which says in essence that the proportionality constant is the same for both types of forces; dissipation without fluctuation is forbidden.NETDSM (); STD2 () (The violation of this principle is a significant weakness of the NEMD method and of other molecular dynamics algorithmsTodd07 (); Hoover08 () that use solely a deterministic thermostat.) The fluctuation-dissipation theorem guarantees the preservation of the non-equilibrium probability in phase space, and hence that time averages equal phase space averages.

The fluctuation-dissipation parameter is most simply taken to be a constant, . However, the analysis remains formally exact if it varies with particle, . As a large class of non-equilibrium systems consist of boundary driven flow, including the present shear flow, it is useful to take this to be dependent on position, . Since this parameter represents the strength of the interaction between the reservoir and the sub-system, it is sensible to make this parameter large near the boundary between the two, and small in the interior of the sub-system. In this work, where the reservoir boundaries are located at , one functional form that is explored below is

 θ(r)=[2zLz]2nθ,n=0,1,2,… (2.6)

For this is constant and the influence of the fluctuation-dissipation parameter is uniform throughout the sub-system. As is increased, its direct influence is increasingly confined to the boundary region.

A second form that is explored below is

 θ(r)={0,|z|≤Lz/2−Rcut,θ,|z|>Lz/2−Rcut. (2.7)

This confines the direct effect of the reservoir to slabs of width next to each boundary. In this work was chosen to be the same as the cut-off to the pair potential, but this is not essential. This form for the fluctuation-dissipation parameter will be referred to as the slab form, or the slab thermostat. Periodic boundary conditions are used for . (In the case of shear flow, these are Lees-Edwards sliding brick boundary conditions.)Lees72 ()

The slab form of position dependence is similar to the simulation geometry used by Khare et al.Khare97 () and by Hoover et al.,Hoover08 () namely the central region evolves adiabatically via Hamilton’s equations of motion. However it does differ in the boundary regions, where in the present case the above equations of motion are used, which means that the particles in the slab undergo reservoir-induced shear flow.

#### ii.1.1 General Constraint

In non-equilibrium thermodynamics, a quantity that often occurs is the adiabatic rate of change of the static part of the reservoir entropy, . The autocorrelation function of this is essentially the Green-Kubo relation for the hydrodynamic transport coefficients.NETDSM () This quantity represents the adiabatic relaxation of relevant macroscopic fluctuations. The adiabatic dynamics can differ significantly when the system is perturbed by the reservoir. In general the reservoir influence affects the dynamic order that this represents, with the result that the transport coefficients can be expected to depend on its value. This expectation has been confirmed for the thermal conductivity AttardIX () and will be demonstrated for the shear viscosity below.

As is discussed in §8.3 of Ref. [NETDSM, ], a general way of circumventing the problem is to impose a constraint on the reservoir force so that it is orthogonal to the gradient of the adiabatic flux,

 ¯¯¯¯¯R⋅∇p˙S0r,st(Γ,t)=0. (2.8)

If is the component of the unconstrained dissipative force given above, then one can introduce a Lagrange multiplier so that the constrained dissipative force is

 ¯¯¯¯Riα=¯¯¯¯Rfiα+μθi∂piα˙S0r,st. (2.9)

The Lagrange multiplier is given by

 μ(Γ,t)=−∑i,α¯¯¯¯Rfiα∂piα˙S0r,st∑i,αθi(∂piα˙S0r,st)2. (2.10)

Three points can be made about this procedure. First, there is no need to add a similar constraint to the stochastic part of the reservoir force because it averages to zero, as has been confirmed numerically. Second, this constraint effectively projects the dissipative force onto a hypersurface of dimension . Since this is negligibly different to the dimensions of the unconstrained dissipative force, imposing it should have negligible effect on the preservation of the non-equilibrium phase space probability density by the equations of motion. And third, if any particle is in a region such that , then the unconstrained , constrained , and stochastic reservoir forces are also zero, which means that the particle moves adiabatically.

#### ii.1.2 Equipartition Theorem

The equipartition theorem, which is formally exact for a non-equilibrium system, is most generally written asNETDSM (); STD2 ()

 ⟨∇∇Sr(Γ,t)⟩=−k−1B⟨[∇Sr(Γ,t)][∇Sr(Γ,t)]⟩. (2.11)

This holds for the position components if, and only if, there is a confining potential that makes the density vanish at and beyond the boundary of the sub-system. If one restricts this to momentum components, takes the trace, and neglects the fluctuations in the consequent extensive quantity, then on the likely points of phase space this becomes

 ∇2pSr(Γ,t)=−k−1B[∇pSr(Γ,t)]⋅[∇pSr(Γ,t)]. (2.12)

For the canonical equilibrium system, . Inserting this yields the equilibrium equipartition theorem, .

As mentioned above, the reservoir entropy of a non-equilibrium system is composed of static and dynamic parts, . By definition, since the static part corresponds to an instantaneous equilibrium system it must have even parity with regard to velocity reversal. Since this is arguably the dominant contribution, one need only retain the odd projection of the dynamic part of the reservoir entropy, .NETDSM (); STD2 () Since the dynamic part vanishes in the equilibrium case, say , it must be at least linear in the non-equilibrium parameter, , . (The non-equilibrium parameter can characterize the strength of an applied thermodynamic gradient, or of an external time-varying potential, as generic examples.) If one restricts attention to dyadic elements of the same parity, say to be definite, then the left hand side of the generalized equipartition theorem in the non-equilibrium case is

 ⟨∇p∇pSr(Γ,t)⟩ne=⟨∇p∇pSr,st(Γ,t)⟩st+O(γ2), (2.13)

and the right hand side is

 −k−1B⟨[∇pSr(Γ,t)][∇pSr(Γ,t)]⟩ne = −k−1B⟨[∇pSr,st(Γ,t)][∇pSr,st(Γ,t)]⟩st −k−1B⟨[∇pSoddr,dyn(Γ,t)][∇pSoddr,dyn(Γ,t)]⟩st+O(γ2) = −k−1B⟨[∇pSr,st(Γ,t)][∇pSr,st(Γ,t)]⟩st+O(γ2).

A similar result holds for elements. Taking the trace of this (which makes it extensive), and neglecting fluctuations (because for an extensive quantity they are relatively negligible), gives on the likely points in phase space

 ∇2pSr,st(Γ,t) = −k−1B[∇pSr,st(Γ,t)]⋅[∇pSr,st(Γ,t)]+O(γ2).

This is the non-equilibrium analogue of the usual equilibrium equipartition theorem; in both cases the left hand side is proportional to the reservoir temperature.

It is important to emphasize that this result is only correct to zeroth and linear order in the non-equilibrium parameter. The point is that many NEMD simulation algorithmsTodd07 (); Hoover08 () invoke the equipartition theorem to obtain the temperature of the configuration from the peculiar kinetic energy, which is then used in the thermostat. Such an argument is based on a locally co-moving Lagrangian reference frame. However, since it neglects the contribution from the velocity gradient, this expression cannot be used for strongly non-equilibrium systems because it is not exact beyond leading order. This is one of the reasons for the uncertainty in existing simulations for the behavior of the shear viscosity at high shear rates.

Beyond the equipartition theorem, Rugh has given a rather general expression for the temperature as a configurational average for an isolated, equilibrium system. Rugh97 (); TDSM () This has also been used in NEMD simulations,Delhommelle03 () even though there is no justification for it as a non-equilibrium average. In fact, since the usual equilibrium equipartition theorem is a special case of Rugh’s expression, the result given above demonstrates that Rugh’s expression cannot be valid beyond leading order in the non-equilibrium parameter. (A critical discussion of the application of a Rugh temperature to non-equilibrium systems is also given by Hoover and Hoover.)Hoover09 ()

Beyond the above criticism (that neither the equipartition theorem nor Rugh’s expression are applicable beyond linear order in the non-equilibrium parameter), it should also be pointed out that the left hand side gives the reservoir temperature, not the sub-system temperature. In an equilibrium system these are equal, of course, and one can regard the equipartition theorem as giving the sub-system temperature at each phase space point. Because it is an equilibrium system, one can further demand that this sub-system configurational temperature be uniform and equal to the reservoir temperature, which makes it useful as a thermostat control function. In a non-equilibrium system neither property is true. One cannot assume that the sub-system temperature is equal to the reservoir temperature, and one cannot assume that the sub-system temperature is uniform throughout the sub-system. Both assumptions become increasingly dubious as the non-equilibrium parameter is increased, and they create doubts about any thermostatted equations of motion in the non-linear regime. (The present SDMD equations of motion invoke only the reservoir temperature, and have no need for a sub-system configurational temperature, and make no assumptions regarding its homogeneity.)

### ii.2 Shear Flow

#### ii.2.1 Isolated System

The velocity is the thermodynamic conjugate of the momentum (see §9.6 of Ref. NETDSM, ),

 ∂σ(r)∂p(r)=−1T(r)v(r). (2.16)

Here is the (first) entropy density of an isolated system, , is the momentum density, is the particle mass, is the number density, is the velocity, and is the temperature.

Based on this, for shear flow one formulates macrostates as specified values of the zeroth and first sub-system momentum moments, to which the zeroth and first sub-system velocities (see §II.2.2 below for their interpretation) are thermodynamically conjugate,NETDSM ()

 vs0≡−T∂Ss(P0,P––––1)∂P0 and v––––s1≡−T∂Ss(P0,P––––1)∂P––––1. (2.17)

The other state variables are not shown explicitly. The zeroth and first momentum moments are

 P0≡∫Vdrp(r) and P––––1≡∫Vdrrp(r). (2.18)

From momentum conservation, the rate of change of the momentum density is just the negative divergence of the flux. Hence the adiabatic rate of change of the first momentum moment tensor is

 ˙P––––01 = ∫Vdrr˙p0(r) (2.19) = −∫Vdrr∇⋅J––––(r) = = VJ––––.

This assumes that the momentum flux is uniform over the volume and vanishes on the boundary.

In general the total momentum flux tensor consists of the pressure tensor and the convective momentum flux, (see §5.1.4 of Ref. NETDSM, ). The pressure tensor comprises the thermodynamic pressure and the viscous pressure tensor, , with the latter also known as the diffusive momentum flux, . One sees therefore that the adiabatic rate of change of the first momentum moment with the convective momentum flux removed is proportional to the pressure tensor.

In the phase space of the sub-system, the zeroth and first momentum moments are

 P0(Γ)=N∑i=1pi, and P––––1(Γ)=N∑i=1qipi. (2.20)

The zeroth moment is a conserved variable, and hence its adiabatic rate of change is zero. The adiabatic rate of change of the first momentum moment is

 ˙P––––01(Γ) = N∑i=1[˙q0i⋅∇qiP––––1(Γ)+˙p0i⋅∇piP––––1(Γ)] (2.21) = N∑i=1[˙q0ipi+˙p0iqi] = N∑i=1[pipim+Fiqi].

Defining the peculiar momentum as , this may be written as

 ˙P01––––––(Γ)=N∑i=1[~pi~pim+Fiqi]+mN∑i=1v(qi)v(qi). (2.22)

Removing the convective momentum flux as discussed above, the pressure tensor is therefore

 P––––(Γ) = 1V˙P01––––––(Γ)−mVN∑i=1v(qi)v(qi) (2.23) = 1VN∑i=1[~pi~pim+Fiqi] = 1VN∑i=1~pi~pim−1VN∑i

The final equality holds for a pair-wise additive potential and it shows explicitly the symmetry of the pressure tensor. The Irving-Kirkwood stress tensor is the negative of this.

The linear constitutive relation from hydrodynamics relates the most likely value of the traceless part of the viscous pressure tensor to the traceless symmetric part of the velocity gradient tensor,

 ¯¯¯¯Π––––∗=−2η[∇v(r)]∗,sym≡−2ηv––∗,sym1, (2.24)

where is the shear viscosity. In component form this is

 ¯¯¯¯Π∗αβ=−η[∂vα(r)∂rβ+∂vβ(r)∂rα−23δαβ∇⋅v(r)]. (2.25)

This result holds in a most likely sense, going forward in time.NETDSM ()

#### ii.2.2 Reservoirs

Consider two reservoirs beyond boundaries located at , moving with velocities . The zeroth and first velocities exerted by the reservoirs on the sub-system areNETDSM ()

 vrx0 ≡ 12[vrx++vrx−] and γ≡vrx1 ≡ 1Lz[vrx+−vrx−]. (2.26)

The zeroth velocity is the average or mid velocity, and the first velocity is the velocity gradient or applied shear rate . It is emphasized that the applied shear rate (as well as the applied temperature that enters below) is a property of the reservoirs. As such it is independent of the state of the sub-system, which means that it has the same value if the resultant sub-system motion is steady shear flow or if it undergoes unsteady or turbulent flow. This observation means that using or in the present equations of motion does not assume a particular sub-system profile, nor does it introduce any unintended bias into those equations of motion.

The phase space probability distribution is essentially the exponential of the reservoir entropy, . Since the reservoir velocities are equal and opposite, and . In this case, since energy and momentum are exchangeable with the reservoir, and since velocity is conjugate to momentum, the static part of the reservoir entropy is

 Sr,st(Γ) = −H(Γ)T+γTPzx1(Γ). (2.27)

The first term on the right-hand side arises from the exchange of energy with the reservoirs, and the second term from the exchange of momentum. As mentioned, both the temperature and the shear rate that appear here are properties of the reservoir that are independent of the state of the sub-system. Likewise the product is dictated by the spatial arrangement of the reservoirs and one does not have the thermodynamic freedom to replace this by some local quantity that depends on the state of the sub-system.

With this expression for the static part of the reservoir entropy for shear flow, the components of the unconstrained dissipative force are

 ¯¯¯¯Rfiα(Γ) = θiα|Δt|2kB∂Sr,st(Γ)∂piα+θiα|Δt|2kB(ˆt−1)¯¯¯¯¯¯¯¯¯¯S′piα (2.28) = −θiα|Δt|2kB[piαTmi−γTziδαx] −θiα|Δt|2kB(ˆt−1)[¯¯¯piαTmi−γT¯¯¯ziδαx] = −θiα|Δt|2mikBT[piα−miγziδαx].

One can recognize the bracketed term as the peculiar momentum, which is on average zero (to ). The constant term in the first and second equalities vanishes for a forward time step, since . The final result is of the same form as a hydrodynamic drag force in the Langevin equation, but using the local fluctuating velocity of the atom.Kremer90 (); Shang17 () Unlike the Langevin equation and its peculiar extension, the present result has a rigorous thermodynamic and statistical mechanical justification.

#### ii.2.3 Constraint for Shear Flow

As was discussed above it is useful to add a constraint to the dissipative reservoir force to make it orthogonal to the gradient of the adiabatic rate of change of the static part of the reservoir entropy. In the case of shear flow, the latter is proportional to . Hence the constraint to be applied is

 ¯¯¯¯¯R(Γ)⋅∇p˙P0zx1(Γ)=0. (2.29)

Using the expression given above, one has explicitly

 ∂˙P0zx1(Γ)∂piα=m−1i[pxiδzα+pziδxα]. (2.30)

With the unconstrained dissipative force given above, one can introduce a Lagrange multiplier so that the constrained dissipative force is

 ¯¯¯¯Riα=¯¯¯¯Rfiα+μθiαm−1i[pxiδzα+pziδxα]. (2.31)

The Lagrange multiplier is given by

 μ(Γ,t)=−∑i,α¯¯¯¯Rfiαm−1i[pxiδzα+pziδxα]∑i,αθiαm−2i[pxiδzα+pziδxα]2. (2.32)

Obviously both the free and the constrained dissipative force, as well as the stochastic force , vanish for molecules in regions where .

#### ii.2.4 Equipartition Theorem for Shear Flow

For shear flow, the equipartition theorem, Eq. (II.1.2), with the above result for the static part of the reservoir entropy, Eq. (2.27), may be rearranged as

 3NmkBT=∑iα[piα−mγziδxα]2+O(γ2). (2.33)

The terms that appear explicitly here may be recognized as being the same as for the equilibrium equipartition theorem, but for the peculiar kinetic energy rather than the total kinetic energy. One can use this to define a local kinetic temperature, , but this is not equal to the actual sub-system temperature except at low shear rates.

### ii.3 Viscous Heating

In shear flow, viscous heating raises the temperature of the sub-system until it is balanced by heat conduction to the boundary. The rate of doing work is the velocity times the force, and the latter is the rate of change of momentum. The rate of work performed on a slab at is

 ˙W = −A[¯¯¯¯Πxz¯¯¯vx(dz/2)−¯¯¯¯Πxz¯¯¯vx(−dz/2)] (2.34) = Adzηγ2.

The negative sign is because is the momentum leaving the slab through the upper face. With this, the rate of change of temperature is

 ˙T(z) = c−1p[˙w(z)−1A∇⋅¯¯¯J0E(z)] (2.35) = c−1p[ηγ2+λ∇2T(z)],

where is the heat capacity per unit volume and is the thermal conductivity. In a steady state system the temperature is constant, , and so

 T(z) = −ηγ22λz2+const. (2.36) = T+ηγ22λ[L2z4−z2].

The second equality holds if the temperature at the boundary equals that of the reservoir. This expression is valid in regions where the heat being extracted by a thermostat is negligible (ie. regions of adiabatic evolution). Similar expressions have been obtained by Khare et al.Khare97 () and by Hoover et al.Hoover08 ()

This variation in temperature is obviously significant for large shear rates and large systems. This is one source of non-Newtonian behavior, and, depending on the sign of (both signs occur in the results below), it will contribute to thickening or thinning with increasing shear rate. This also causes the viscosity to depend on the system size at high shear rates, which is extremely unusual in thermodynamics, since, at least in the equilibrium case, in the thermodynamic limit all parameters are either intensive (ie. independent of system size) or else extensive (ie. linearly proportional to the system size).

This temperature dependence and non-intensive behavior poses conceptual challenges in how to apply results for the shear viscosity. It is emphasized that this is a real physical effect, since in steady state shear heat can only be dissipated via the boundaries.

One pragmatic solution is to perform simulations in the regime where the temperature increase is relatively negligible, , or

 γ≪√8λTηL2z or |vr,+−vr,−|≪√8λTη. (2.37)

If this holds, we may assume that the sub-system temperature is the same as the reservoir temperature. In this case the shear viscosity can be expected to be independent of system size, . In this regime the temperature of the sub-system is the same as the temperature of the reservoir, and there is no ambiguity as to which temperature is referred to.

In the real world or laboratory situation, very high shear rates only occur over a limited spatial extent, such as in a shock wave, or in the vicinity of an injection aperture, or around a sharp protuberance. In addition, very high shear is often of a transient nature, and the amount of viscous heat produced may be insignificant compared to the heat capacity of the relevant volume. In these cases one can assume that the temperature in the shearing region is constant, and one can use the results for the shear viscosity that are obtained in the limit where it is independent of the system size, and where the temperature of the sub-system is approximately uniform and equal to that of the reservoir.

### ii.4 Pressure

The normal component of the thermodynamic pressure is now analyzed. With the applied shear rate being , the derivative can be made at constant shear rate or else at constant reservoir velocities . The second of these is arguably the case with greatest physical relevance. One can define the -pressure as the derivative of the total unconstrained entropy,

 pz ≡ TA∂Stot(N,V,T,γ)∂Lz (2.38) = kBTAZtot∂∂Lz∫dΓeSr(Γ)/kB = kBTAZtot{3NLzZtot+∫dΓeSr(Γ)/kB∂Sr(Γ)kB∂Lz} = 3NkBTV+TA⟨∂Sr(Γ)∂Lz⟩ne.

The first term is that of an ideal gas, and it follows from the usual re-scaling of the sub-system. As mentioned above the reservoir entropy is the sum of static and dynamic parts, with the static part being given by Eq. (2.27), . Writing , and supposing that the potential energy is the sum of pair potentials, one has

 ∂H(Γ)∂Lz = ∑i

Also, with the first momentum moment being , a similar scaling yields

 ∂Pzx1(Γ)∂Lz={L−1zPzx1(Γ),γ=const.,0,vr,±=const. (2.40)

Hence one has (for the case )

 pz = 3NkBTV−1V⟨∑i

The first two terms may be recognized as essentially the non-equilibrium average of the usual equilibrium virial expression for the pressure. Since , the third term for the case is

 γV⟨Pzx1(Γ)⟩ne = γVmγN∑i=1⟨z2i⟩ne+O(γ4) (2.42) = mnγ2L2z12+O(γ4).

The final approximation assumes a uniform density profile. One sees that there is a repulsive contribution to the pressure that goes as . This contrasts with a uniform equilibrium system where the pressure is intensive, which is to say that it is independent of . This situation is similar to that of an inhomogeneous equilibrium fluid (since a shearing system is also inhomogeneous, Eq. (2.36)), such as a slit pore, where the normal component of the pressure also depends on separation . In a sense it is worse in the present case in that the pressure actually diverges quadratically with sub-system width.

As mentioned this derivative at constant shear rate is different to the derivative at constant reservoir velocity . In the latter case is constant and the derivative of the first momentum moment vanishes. In this arguably more physical case, the normal component of the pressure does not have such a divergent non-extensive contribution. It likely still has an dependence due to the inhomogeneous nature of the shearing system.

### ii.5 SDMD Simulation Details

A Lennard-Jones fluid was simulated, with cut and shifted pair potential

 (2.43)

Here , but in fact this is immaterial. All of the following results are made dimensionless using the well-depth energy , the atomic diameter , and the mass . These give the unit of time as , which is 2.2 ps for argon.

The cut-off radius was . It can be shown that the consequently neglected tail contribution to the shear viscosity is , which is negligible. This was confirmed numerically by performing one simulation with , where it was found that the shear viscosity changed by 0.09%, which is insignificant compared to the 0.4% statistical standard error for those runs.

A rectangular simulation cell was used, usually with . In all cases 1000 atoms were used. Periodic boundary conditions were applied, with Lees-Edwards sliding brick conditions imposed in the -direction.Lees72 ()

A small cell spatially based neighbor table was used, with cells cubic and of edge length 0.7.AttardV (); AttardIX () This typically gave a neighborhood volume 2.6 times the cut-off sphere, which is much better than typical large cell neighbor tables achieve. The neighbor table was modified for the Lees-Edwards conditions by adding a separate neighbor table of width beyond the upper and lower boundaries and maintaining a record of the images of the atoms in this region. This extension to the neighbor table was fixed to the central neighbor table, and the image atoms moved through it in accord with the Lees-Edwards conditions. The height of the central cell was necessarily an integer multiple of the small cell edge length, but not necessarily so for the extended system . The neighbor table performed remarkably well, consistent with previous experience where the simulations scaled linearly with the number of atoms. AttardV (); AttardIX ()

The Verlet leap frog algorithm was used.Shang17 (); Gao16 () This is second order in the time step, and is about an order of magnitude more efficient than simple time stepping. The time step was , which was halved in some tests. Higher order methods were not explored.Shang17 ()

Usually simulations were performed for time steps, with high shear rates requiring fewer time steps. To obtain the temperature derivative of the viscosity time steps were used in total at each state point. Each case was usually started from a nearby shear rate, after velocity re-scaling and position pinching, and a brief equilibration period. The neighbor table was updated every time step. Values for averaging were accumulated once every 20 time steps. In hindsight, it would be more efficient to accumulate these more frequently. Each simulation was broken into 50 blocks, the standard deviation of which were estimated from the fluctuations between them. The statistical error quoted everywhere below is one standard error on the mean, which is this standard deviation divided by . It gives the 68% confidence interval.

## Iii Results

Figure 1 shows the viscosity as a function of the fluctuation-dissipation parameter. The state point was , , and for these simulations, , , and . The shear rate is fixed at , which corresponds to 45 GHz for argon. Based on results below at this state point, this can be considered to be in the linear, low shear rate regime. The limit for negligible temperature change, Eq. (2.37), is , using AttardV () and .

First focussing on the uniform application of the fluctuation-dissipation parameter, , and the unconstrained equations of motion (empty triangles), one sees that the simulated viscosity is strongly affected by the strength of the parameter, with the viscosity decreasing as the parameter is increased. It will be recalled that the stochastic, dissipative equations of motion are correct to linear order in , which means that they are unreliable for large values. Enforcing the constraint, Eq. (2.29) (filled triangles), largely removes this sensitivity.

The influence of on the viscosity is also reduced when it is focussed in the boundary region, (squares), (diamonds), and confined to the boundary slab (circles). In these cases the viscosity (and the temperature) are measured in the central half of the system, here and below. In the slab case the evolution in this central region is strictly adiabatic. One can see that at , the unconstrained results (open symbols) increase toward the constrained, uniform result as the dissipative and stochastic forces are more narrowly confined to the boundary. There is relatively good agreement between the different boundary methods. The slight disagreement between the constrained and unconstrained boundary slab results at indicates that the structure and motion in the boundary slab region can still influence the viscosity measured in the central region.

If the fluctuation-dissipation parameter is too small to remove the viscous heat, then the sub-system kinetic temperature will be larger than the reservoir temperature, particularly for a boundary-confined parameter. This will effect the viscosity even in the constrained case. For example, for the constrained case, at , the central viscosity was and the central kinetic temperature was . The temperature derivative at this state point was measured to be (measured at with and ). Using this to correct the measured viscosity to gives , which is closer to the values measured at , . The correction is not exact because of additional effects from the induced inhomogeneity in both temperature and density.

This particular state point, , , was chosen to make quantitative comparison with the viscosity simulated with four different methods by Hess.Hess02 () Hess obtained from the transverse current method and , from the pressure fluctuations (Green-Kubo) method , from a periodic external perturbation method , and from the NEMD (SLLOD, ) method . These are in agreement with each other but they are significantly smaller than the present best estimate of the shear viscosity, 0.54–0.55. However, the four types of simulation performed by Hess are not independent of each other as they were all deterministic molecular dynamics using a uniform Berendsen thermostat with the same coupling time of 20 (all methods).Hess02 () Since this is presumably large enough to hold the kinetic temperature at the specified reservoir value, the most direct comparison with the present results is the uniform, , unconstrained case, with being similarly large enough. This case gives , which is comparable to Hess’ estimates. (As mentioned in §II.5 above, the cut-off radius, here, in Ref. [Hess02, ], has negligible effect on the viscosity.)

Figure 2 shows the local kinetic temperature, Eq. (2.33). It can be seen that the uniform fluctuation-dissipation parameter, , , is relatively good at maintaining a uniform kinetic temperature, with average value . The difference from the reservoir temperature appears to correspond to the neglected corrections to the kinetic temperature. For the boundary slab case, the kinetic temperature profile is parabolic in the central adiabatic region, with the curvature given by Eq. (2.36), using the thermal conductivity (interpolated from Ref. [AttardV, ]) and the present viscosity . The kinetic temperature profile (circles) has slightly larger curvature than the macroscopic temperature profile (curve), which may either be interpreted as a small change in the value of the shear viscosity or the thermal conductivity, possibly because of the inhomogeneities, or else as weak higher order contributions to the equipartition theorem.

Figure 3 shows the density profile for the uniform and slab parameter cases. In the uniform case, the density depletion with increasing distance from the center appears to correlate with the increase in the mean speed and hence the kinetic energy . The explanation is that it is entropically favorable to transfer peripheral particles to the center where they have lower speed because this releases energy to the reservoir and hence the entropy associated with the phase space point is increased. For the uniform parameter case, it was found that the peak of the density profile increased with increasing at constant . This is consistent with the slab result, where the density is enhanced in the slab region compared to the adiabatic region. Undoubtedly these density inhomogeneities contribute to the sensitivity of the shear viscosity to the value of the fluctuation-dissipation parameter. Given the 1–2% difference in temperature and density profiles between the uniform and slab application of the fluctuation-dissipation parameter, it is unrealistic to expect the two methods to produce viscosities in better agreement than this. Ultimately the two formulations model different physical situations, with the uniform case applicable to the small size asymptote, and the boundary slab applicable to the size-dependent adiabatic case.

Figure 4 shows the viscosity as a function of shear rate at the same state point as above, , . The constrained equations of motion were used with uniform (, triangles) and boundary slab (circles) application of the fluctuation-dissipation forces. Except at the highest shear rates, the same value of the fluctuation-dissipation parameter was used. At low shear rates the uniform and the slab results asymptote to the same constant viscosity, . It can be seen that for the uniform case, there is a noticeable shear thinning for shear rates . Recall that for this state point the limit for negligible temperature change, Eq. (2.37), is .

For the boundary slab case, the shear viscosity was obtained by dividing the average component of the diffusive momentum tensor in the central half of the sub-system by the gradient in the velocity, , actually in that region. At this and the higher densities below, the latter was almost always higher than the applied shear rate, increasingly so as the applied shear rate was increased.

For the boundary slab case, the kinetic temperature in the adiabatic central region increased significantly at high shear rates, and a correction was used to facilitate comparison with the uniform fluctuation-dissipation parameter case. The temperature corrected viscosity is , where the temperature derivative is obtained by finite difference of two simulations at . In this case was obtained at with and . Also, is the average viscosity and is the average kinetic temperature in the central adiabatic half.

This temperature correction is problematic at high shear rates for three reasons. First, at high shear rates the kinetic temperature is not the actual temperature, Eq. (2.33). Second, no correction is made for the density difference. And third, the linear correction is unreliable for large temperature shifts.

Despite these concerns, it is clear that the temperature difference accounts for much of the difference between the viscosity measured in the uniform and in the slab cases at low to moderate shear rates, . At the highest shear rate , the actual viscosity in the slab case is larger than that in the low shear rate limit, whereas the temperature-corrected viscosity is smaller than the limiting value. One could describe the raw data as shear thickening, and the corrected data as shear thinning. In this case it is clear that the shear thickening comes entirely from the temperature increase due to viscous heating in the adiabatic region, since at this state point .

Figures 5 and 6 shows the viscosity as a function of shear rate at and at and , respectively. The limiting shear rate for negligible temperature increase, Eq. (2.37), is for (using AttardV () and ), and for (using AttardV () and ).

Again both the uniform and the slab fluctuation-dissipation methods agree on the low shear rate asymptotes. The low shear rate viscosity evidently increases with increasing density. Further both methods, at least initially, predict shear thinning for both densities, although the quantitative agreement between the two models is limited at higher shear rates.

At , the raw boundary slab results show shear thickening at higher shear rates due to the increasing temperature in the adiabatic central region and the positive value of the temperature derivative of the viscosity. Correcting for this temperature effect gives shear thinning behavior.

Interestingly enough, at the higher density the temperature derivative of the viscosity is negative, and both the raw and the corrected results show shear thinning. The temperature corrected boundary slab results again lie closer to the uniform fluctuation-dissipation parameter results, which confirms that much of the difference between the two can be attributed to the higher temperature in the adiabatic region.

For most of the results at , increasing the value of in the uniform constrained case increased the value of at a given . This is the opposite behavior to that observed in the above results at the lower density. Possibly the sign of this effect depends on the sign of the temperature derivative of the viscosity, since increasing generally decreases the kinetic temperature.

All of the results presented above and below are for a time step of . It should be mentioned that halving this to at and changes the various averages by 0.1–0.5%, which is about an order of magnitude larger than the typical statistical error. This systematic error is less at lower shear rates.

The quantitative effects of finite size were estimated at , , and , in the uniform case (1000 particles). At , changing from 11.9 to 15.4 decreases by 0.2% and decreases by 1%. Changing from 5 to 10 at changes from 2.39 to 2.20, and decreases by 5%.

Density profiles for the uniform case are shown in Fig. 7. The density profile is of the form , with the central density being independent of system size and increasing with increasing fluctuation-dissipation parameter. The corresponding temperature profiles are uniform, show a similar independence from system size, and decreases with increasing fluctuation-dissipation parameter value. More broadly, the magnitude of the density inhomogeneity increases with shear rate and with fluctuation-dissipation parameter.

The kinetic temperature profiles shown in Fig. 8 confirm the uniformity for the case, and the increasing efficacy of the fluctuation-dissipation parameter with increasing magnitude. For the boundary slab case, one again sees that the continuum expression Eq. (2.36) is applicable in the central adiabatic region. The small discrepancy between this expression and the measured values can either be taken to indicate small local changes in or as a result of the inhomogeneity in or , or else it indicates that the difference between the kinetic temperature and the actual temperature, the neglected higher order contributions in Eq. (2.36), are important at this high shear rate.

The density profiles shown in Fig. 9 are qualitatively similar to those shown for the lower density case, Fig. 3. The sharp density depletion at the periodic boundary in the uniform fluctuation-dissipation parameter case, , induces an adjacent clear local maximum or density oscillation at this high density. The behavior of for the boundary slab case is qualitatively different, with a minimum at . It can be argued that the density enhancement in the slab region is consistent with the uniform parameter case in that a large fluctuation-dissipation parameter apparently favors a higher density.

At , , for the boundary slab case with , in the central adiabatic region the simulations show that , , and . (The change upon halving the time step to in this case is smaller than the statistical error.) Evidently then . This ordering agrees with Hoover et al.,Hoover08 () and it appears significant compared to any possible systematic error. Also, for the pressure tensor in the same region, , , and . Hence . This ordering agrees with Hoover et al.Hoover08 () Halving the time step to changes these to , , and .

At the same state and shear point for the uniform case , , the simulations show that for the whole system , , and