The viscoelastic behavior of sheared fluids is calculated by Non-Equilibrium Molecular Dynamics (NEMD) simulation, and complementary analytic solutions of a time-dependent extension of Eyring’s model (EM) for shear thinning are derived. It is argued that an “incremental viscosity,” *η*_{i}, or IV which is the derivative of the steady state stress with respect to the shear rate is a better measure of the physical state of the system than the conventional definition of the shear rate dependent viscosity (i.e., the shear stress divided by the strain rate). The stress relaxation function, *C*_{i}(*t*), associated with *η*_{i} is consistent with Boltzmann’s superposition principle and is computed by NEMD and the EM. The IV of the Eyring model is shown to be a special case of the Carreau formula for shear thinning. An analytic solution for the transient time correlation function for the EM is derived. An extension of the EM to allow for significant local shear stress fluctuations on a molecular level, represented by a gaussian distribution, is shown to have the same analytic form as the original EM but with the EM stress replaced by its time and spatial average. Even at high shear rates and on small scales, the probability distribution function is almost gaussian (apart from in the wings) with the peak shifted by the shear. The Eyring formula approximately satisfies the Fluctuation Theorem, which may in part explain its success in representing the shear thinning curves of a wide range of different types of chemical systems.

## I. INTRODUCTION

Theoretical treatments of the shear viscosity, *η*, of liquids started with empirical expressions for its temperature and pressure dependence. The effect of shear is to distort the molecular arrangement of the molecules which induces an opposing shear stress that is non-linear in the shear rate (SR) when the shear rate has large values. This gives rise to shear thinning, which has also been represented by semi-empirical analytic expressions. The molecular level factors that govern the shear rate dependent viscosity are still poorly understood,^{1,2} and while the distortion of the local structure of a liquid under shear is now well characterized for small molecules,^{3,4} theoretical prediction of the shear rate dependence of the viscosity^{5} has proved to be more problematic. Phenomenological models of shear thinning are still prevalent in part for this reason.^{6,7}

In 1936, Eyring derived an expression for the shear rate dependence of *η* based on what today might be called a “trap” model,^{8} which is based on assumed processes that occur at the molecular level. This Eyring model (EM) considers a representative molecule to be in a potential well whose escape barrier height is lower in the flow direction owing to the anisotropic structural distortion of the surrounding liquid caused by the shear. The nonlinearity of the stress *vs.* shear rate in the EM leads to the prediction of shear thinning. The temperature and pressure dependence of the Newtonian viscosity part can be included largely as a separate exercise in the expression for the prefactor. The appealing feature of this model is that the assumed origins of the shear thinning are intuitively reasonable and the parameters could in principal be related to molecular properties, therefore making it potentially molecule specific to a certain extent. The EM for shear thinning has been used widely in fitting to experimental viscosity data, such as in high pressure lubrication,^{9} and even for solids undergoing plastic flow.^{10,11} We note in passing that more recently, other trap models have been derived which are more relevant to slowly evolving or ageing “quenched” systems such as glasses and gels.^{12}

Nonequilibrium statistical mechanics has developed over recent decades to provide useful analytic expressions relating to the shearing of liquids, some aspects of which can be implemented in molecular dynamics simulation. In the low shear rate limit, the Green-Kubo equation,^{13} or the equivalent Helfand-Einstein formula,^{14} can be used to determine the viscosity. Another significant theoretical development describes the stress build-up from an equilibrium state after the application of a shear rate of arbitrary magnitude, which is cast in the form of a *transient time correlation function* (TTCF).^{15–17} In 1993, Evans, Cohen, and Morriss proposed a formula expressing the probability of states for a system driven arbitrarily far out of equilibrium, which is called the Fluctuation Theorem (FT).^{18–20} In fact there are many slight variants of the FT formula which reflect different conditions applied to the system.^{15,16,20,21} Non-Equilibrium Molecular Dynamics (NEMD) simulation is a complementary tool which can be used to compute these functions for model fluids and thereby provide definite predictions for specific chemical systems.

Constitutive relations for shear viscoelasticity in principle provide a complete description of the time dependent rheology of a liquid and have been widely investigated.^{22–24} This work focuses on the steady-state stress relaxation function, *C*(*t*), or SRF which requires a different procedure to compute by NEMD, and has not been investigated far from equilibrium by this technique, as far as we are aware. The steady state shear SRF currently falls outside the scope of the above mentioned non-equilibrium statistical mechanical advances. The SRF is the transient viscoelastic response of a sheared liquid to a suddenly applied additional infinitesimal shear rate and gives directly what might be called an *incremental viscosity*, which could be considered to be a better measure of the system at a given shear rate than the usual definition of the shear-rate dependent viscosity (i.e., the shear stress divided by the shear rate). This is discussed in Sec. II. It is shown how these quantities can be obtained by NEMD and also within the framework of the EM. *Inter alia* analytic expressions are derived for the SRF and TTCF of the Eyring model. In the case of the SRF, insights are gained into why it has proved to be so successful in representing the shear thinning behavior of a diverse range of chemical systems. The extent to which the Eyring model is consistent with the FT is also explored here. It is shown that one of the limitations of the viscoelastic Eyring model derived here is that it does not allow for a time delay in structural evolution after the start of shearing, which could give shear stress overshoot (as is found in a number of experimental systems). A simple modification of the derived equations which has the same consequence as a “fictive” quantity to implement this time-delay is proposed which remedies this problem in a semi-empirical way.

To date, the fields of NEMD and associated non-equilibrium statistical mechanics of fluids and rheological constitutive equations have been advanced largely in parallel and independently. The present study is a preliminary attempt to analyze the links and potential synergy between the two approaches, with the ultimate objective of formulating constitutive equations that are consistent and up to date with the advances and conclusions derived from NEMD and related non-equilibrium statistical mechanics developments.

## II. INCREMENTAL VISCOSITY AND RELATED RELAXATION FUNCTIONS

In this section, the concepts of an incremental viscosity, and related time dependent quantities for an arbitrarily large shear rate, are introduced. Their calculation and investigation using NEMD and within the theoretical framework of the Eyring model are also discussed, particularly in regard to the constraints required by the Fluctuation Theorem (FT) in the latter case.

The viscosity of a fluid in the zero shear rate limit, called the Newtonian viscosity, *η*_{0}, is defined by

where *σ* is the time averaged shear stress at shear rate $\gamma \u0307$ carried out at steady state (the sign convention is such that *σ* > 0 if $\gamma \u0307>0$). In experiment, NEMD, and many theories, the viscosity at a finite shear rate is arbitrarily defined to be $\eta (\gamma \u0307)=\sigma /\gamma \u0307$, which is inconsistent with the definition of Eq. (1) as it is not based on the incremental change in steady state shear stress in response to a small change in the shear rate. A potentially more useful alternative definition of the viscosity at finite shear rates which is consistent with Eq. (1) at any shear rate is what might be called the *incremental viscosity*, *η*_{i}, defined by

which is a more sensitive indicator of the physical state of the liquid at $\gamma \u0307$. Figure 1 gives a schematic diagram illustrating the difference between *η* and *η*_{i}. As shown in Eq. (2), the shear rate dependent viscosity, $\eta (\gamma \u0307)$, is an accumulated average of the incremental viscosities exhibited by the system up to that value of the shear rate and is therefore not a specific measure of the state of the system at that final shear rate, $\gamma \u0307$, unlike the incremental viscosity, $\eta i(\gamma \u0307)$. The incremental viscosity can be obtained from $\eta (\gamma \u0307)$ through

As (usually) $d\eta /d\gamma \u0307<0$, the incremental viscosity is less than the viscosity (i.e., *η*) at a given shear rate, a trend which is accentuated with increasing shear rate, as may be seen directly from the final identity in Eq. (3). The formula is a restatement of the formal definition of *η*_{i}, as $\sigma =\gamma \u0307\eta $, which indicates that the incremental viscosity quantity is definable at any shear rate and is strongly dependent on $d\eta /d\gamma \u0307$.

Some liquids exhibit an increase in shear viscosity after an initial shear thinning stage, an effect known as shear thickening. The viscosity, *η*, has a minimum at this transition from shear thinning to shear thickening.^{25} Equation (3) indicates that when $d\eta /d\gamma \u0307$ changes sign from negative to positive, the incremental viscosity goes from being below to being above the value of *η*. The incremental viscosity definition therefore indicates clearly the shear thinning to thickening transition. The incremental viscosity cannot be zero, *except* in a hypothetical infinite shear rate limit, where *η* would also tend to zero.

So far the analysis has been for systems at steady state. In Sec. II A, the associated time dependent or viscoelastic properties associated with the incremental viscosity definition are considered.

### A. Viscoelasticity

A phenomenological way of representing the linear stress response to a time dependent shear rate is to use the Boltzmann superposition principle,^{26–28} in terms of a time-dependent elastic modulus, *C*, or “shear stress relaxation function” (SRF), which has units of stress,

Equation (4) represents the case where the shear stress is on average zero for *t* < 0, i.e., it is an equilibrium system. In this work, the situation where there is a background or reference shear rate (which can be zero), $\gamma \u0307R$ already established at *t* = 0, is considered. Assume that the system is at steady state at $\gamma \u0307=\gamma \u0307R$ with a time average stress of $\sigma \xafR$. At time *t* = 0, the shear rate is changed to $\gamma \u0307(t)=\gamma \u0307R+\gamma \u0307kf(t)$, which can be considered to be a shear rate “kick,” and where the unitless function, *f*(*t*), takes into account the time dependence of the kick shear rate magnitude, $\gamma \u0307k$,

as the $\gamma \u0307R$ shear rate part of the shear stress is already at steady state, by assuming $\gamma \u0307R$ was initiated sufficiently far in the past, *t* < 0, for the transient start-up response to have finished before *t* = 0. Therefore,

The last step in Eq. (6), which involves replacing $\gamma \u0307$ by $\gamma \u0307R$ in *C*, is made because $\gamma \u0307k$ can always be made arbitrarily small so it does not significantly affect the structural and dynamical state of the system and hence the time dependence of *C*. Therefore we have

where *δ*(*t*′) is the Dirac delta function. The function, $\eta i(t,\gamma \u0307R)\u2192\eta i(\gamma \u0307R)$ in the *t* → *∞* limit, as lim_{t→0}*C*(*t*) = 0. The last equation in Eq. (7) is the corresponding formula in discrete time steps.

For *f*(*t*′) = 1, we also have

The last step derives a widely used relationship in rheology,^{29} which is obtained by making the substitution *x* = *t* − *t*′, so *dx* = −*dt*′. The important point is that far from equilibrium, the appropriate stress relaxation function to be used in the superposition formula is the one which gives *η*_{i} and *not η*, which underscores the importance of the incremental viscosity and related functions. To emphasize this point, the shear stress relaxation function, *C*, used in the kick procedure is denoted by *C*_{i}. In the $\gamma \u0307R\u21920$ limit, *C*_{i}(*t*) tends to the shear stress autocorrelation function of an equilibrium liquid.^{16}

Equations (6) and (8) are key to the methodology used here. The applied shear rate $\gamma \u0307$ is the sum of that of the reference system (which can be zero) and a very small increment of the shear rate. The shear rate increment can be made sufficiently small not to affect noticeably the relaxation characteristics of the system compared to its reference state. The shear stress relaxation function so determined is then that of the reference system, which is what we want. The shear-kick (SK) technique is to probe the viscoelastic state of the reference system. Equation (8) is a useful identity which emphasises a related point. It is shown that this same relaxation function can be obtained when the shear kick is applied continuously during the perturbed trajectory. Both of these approaches are used in the NEMD calculations described in Sec. II B.

Equations (4)–(8) can be written in a form suitable for numerical integration by NEMD by replacing the integrations by summations. For example, in Eq. (4),

where Δ*t* is the time step, *t* = *M*Δ*t* is the time duration of the kicked trajectory, and $\gamma \u0307j$ is the shear rate at time *j*Δ*t*. The formula in Eq. (9) can be used for numerical evaluation of computer generated *C*. Also after the last line in Eq. (7),

The NEMD procedure used here to calculate the incremental viscosity and its corresponding stress relaxation function, *C*_{i}(*t*), for a model simple liquid, and by EM, are presented below.

### B. NEMD

As shown above, the incremental viscosity can be calculated by applying an infinitesimally small shear rate “kick,” $\gamma \u0307k$, at an arbitrary time declared to be *t* = 0 to an equilibrium (i.e., unsheared) fluid, or one that has reached steady state at a constant shear rate, $\gamma \u0307R$. The viscoelastic stress response, *σ*_{k}(*t*), associated with the shear kick (SK) can be written in a form involving the associated shear stress relaxation function, *C*_{i}(*t*), of the fluid in the $\gamma \u0307$ state,

where $\gamma \u0307$ can be 0. The incremental viscosity is $\eta i(t\u2192\u221e,\gamma \u0307)$.

Assuming a monatomic fluid in which the molecules interact with a radially symmetric pair interaction, *ϕ*, the virial expression for the shear stress of the system is

where *N* is the number of molecules in the simulation cell of volume, *V*. For particle *i*, the mass is *m*_{i} (all the masses are the same here) and *r*_{α,ij} is the *α* component of the vector $r\u0332ij=r\u0332i\u2212r\u0332j$ between molecules *i* and *j*. Also, $v\alpha ,i$ is the *α* Cartesian component of the peculiar velocity of molecule *i* (i.e., the excess above any advected component).

A series of trajectories was initiated from equilibrium states or states at steady state shear ($\gamma \u0307$) which were well separated in time. The shear rate kick was either for the duration of a single time step or applied at every time step of the generated new trajectory. The function $Ci(t,\gamma \u0307)$ was obtained directly from the kick lasting just one time step, or alternatively by taking the time derivative of the stress response to a shear rate kick added every time step, and using the right-hand side equation in Eq. (11). The many-step and one step applications of *γ*_{k} correspond to the second and third lines of Eq. (7). This procedure can be used because (a) in the $\gamma \u0307k\u21920$ limit, $Ci(t,\gamma \u0307)$ will in practice be independent of $\gamma \u0307k$ and (b) the reference shear rate, $\gamma \u0307R$, is constant during the kick procedure action. We refer to this small kick variant of NEMD as Kick NEMD or KNEMD, which in slightly different formulations has been applied in the past to model by NEMD shear flow starting from an equilibrium state in Refs. 30 and 31, and by Evans to determine the time dependent structure factor.^{32}

The shear rate of the reference trajectory was $\gamma \u0307=dvx/dy$, where *x* is the flow direction and *y* is the shear gradient direction. Two trajectories starting from the same reference state point were carried out repeatedly (typically 2000 times). Let $q\u03323N$ represent the 3*N* position degrees of freedom and $v\u03323N$ the corresponding quantity for the peculiar molecular velocities (the particle masses are the same). The origin of the coordinate system was at the center of the origin cubic MD cell. Let the microstate set of defining parameters be ${q\u03323N,v\u03323N,xdis,\gamma \u0307R}$, where *x*_{dis} is the Lees-Edwards periodic boundary misalignment distance in the flow direction between adjoining cell boundaries and $\gamma \u0307R$ is the steady state shear rate of the reference system. The initial reference state is labeled Γ_{+}. The reference trajectory and the kicked trajectory with shear rate $\gamma \u0307k$ were carried on from this starting state point. Another reference starting state, Γ_{−}, was obtained from Γ_{+} by implementing *y*-reflection mapping (e.g., see Ref. 16, p. 192), which here is ${qy\u2192\u2212qy,vy\u2192\u2212vy,xdis\u2192\u2212xdis,\gamma \u0307R\u2192\u2212\gamma \u0307R}$, to which an additional very small shear rate, $\gamma \u0307k$, was added to implement the kicked trajectory. The shear rates for either a time step or all time steps in the kicked trajectory were $\gamma \u0307R+\gamma \u0307k$ and $\u2212\gamma \u0307R+\gamma \u0307k$ for the “+” and “−” starting states, respectively. The stress response needed in Eq. (11) is

where *k* stands for a kicked trajectory and 0 stands for one of the two conjugate paired reference trajectories. The stress contributions from the two conjugate pair reference states cancel out and the sum of the two was found to add to zero within machine error for all *t*, a deliberate feature of the procedure which made the number of significant figures defining the kick stress response greater than if conjugate pairing had not been implemented. The magnitude of the kick shear rate should be small enough not to significantly change the fluid’s structure and physical properties. The kick induced stress response function, $Ci(t,\gamma \u0307)$, characterizes or probes the viscoelastic response of the reference fluid and therefore the kick must not perturb this state to any important extent. The stress response must be much smaller than the root mean square stress about the mean stress of the reference fluid.

NEMD simulations of the Lennard-Jones (LJ) fluid were carried out subject to imposed shear flow using velocity rescaling thermostatted SLLOD equations of motion^{33,34} to obtain *η*_{i} and *C*_{i} for that system. The quantities quoted in this work are expressed in the usual LJ pair potential units of energy and distance and the mass, *m*, of the molecule. For compactness of notation, $\gamma \u0307R$ is replaced by $\gamma \u0307$ for the rest of this work as the difference between $\gamma \u0307R$ and the total shear rate during a kick trajectory can be made arbitrarily small.

Figure 2 shows the shear stress relaxation function, *C*_{i}(*t*), associated with *η*_{i} for a LJ liquid at the widely used state point, *ρ* = 0.8442 and *T* = 0.722, which is near the triple point.^{35} Two types of shear history were used to compute $Ci(t,\gamma \u0307)$, (a) by applying a shear rate “pulse,” $\gamma \u0307k=10\u22123$ for 1 time step at *t* = 0 and (b) for $\gamma \u0307k$ applied at every time step for *t* > 0 and using the second equation in Eq. (11) to obtain the SRF. As may be seen in the figure, the two methods give statistically indistinguishable curves for $Ci(t,\gamma \u0307)$ at each value of $\gamma \u0307$. Figure 2 also has in it the time correlation function, *C*_{i}(*t*) = (*V β*)⟨*σ*(0)*σ*(*t*)⟩ (where *V* is the volume of the system, *β* = 1/*k*_{B}*T*, *k*_{B} is the Boltzmann’s constant, and *T* is the temperature), for the *equilibrium* ($i.e.,\u2009\gamma \u0307=0$) reference state,^{36} which is seen to agree very well with that generated by the two kick routes when $\gamma \u0307=0$. The trend is for the *C*_{i}(*t*) to decay more rapidly at long times with increasing the reference fluid shear rate, which is consistent with the observation that the viscosity decreases with increasing shear rate. The value of *C*_{i}(0) increases to a small extent with increasing $\gamma \u0307$. Note also that for $\gamma \u0307=1.5$ and 2.0, *C*_{i}(*t*) is negative at intermediate times, which is not inconsistent with the fact that the stress autocorrelation functions at these shear rates are oscillatory. The figure also shows (green curve) a least squares fit of the correlation function to the Rounded Stretched Exponential (RSE) function defined in Eq. (3) of Ref. 37 and given in the caption of Fig. 2. This analytic form fits well the *C*_{i}(*t*) for $\gamma \u0307=0$. In Ref. 37, the relaxation function is normalised to 1 at *t* = 0, whereas here the *t* = 0 value, *C*_{i}(0), is *G*_{∞} instead.

Figure 3 compares the values of *η* and *η*_{i} as a function of $\gamma \u0307$ for the same system as used for Fig. 2. The figure shows that *η*_{i} decays with increasing shear rate more rapidly than *η*, as expected for a shear thinning liquid. Some of the *η* data in the figure are taken from Ref. 35 for system sizes of mainly 2048 LJ particles. The function,^{38}

where *p*_{0}, *p*_{1}, *p*_{2}, and *p*_{3} are arbitrary variables, can be fitted to the NEMD *η* data quite well, as shown in the figure. The fit parameters are given in the figure caption. The viscosities, $\eta i(\gamma \u0307)$ calculated from Eq. (3) and $\eta (\gamma \u0307)$ from Eq. (14) are also presented in the figure. This analytic IV is seen to be in good agreement with the KNEMD simulation data and within its statistical uncertainty. The agreement seen between the incremental viscosity obtained by (a) just applying the kick shear rate, $\gamma \u0307k$, for the first time step of the segment (labeled “K1” in the figure) and then using Eq. (11) and (b) adding $\gamma \u0307k$ for every time step of the kick (“KA” in the figure) gives strong support for the applicability of Boltzmann’s superposition assumption even when $\gamma \u0307$ is arbitrarily large.

In Sec. II C, exact formulas are derived for the incremental viscosity and related stress relaxation functions for the Eyring model.

### C. Eyring model

The Eyring model (EM) for liquid shear thinning is extended below to encompass the intrinsic viscosity, associated stress relaxation function, and the transient time correlation function (TTCF).

The steady state shear thinning curve of a liquid can be represented by the generic analytic form

where *σ* is the shear stress at a given shear rate, *σ*_{0} is a characteristic shear stress above which the liquid starts to shear thin to a noticeable extent, and *F*(⋯) is a generic function. The term, $\sigma 0$, is a somewhat imprecisely defined quantity and is often obtained by fitting experimental or NEMD simulation shear thinning data to the analytic form given in Eq. (15) and treating *σ*_{0} as a fit parameter. Clearly *F*(*σ*/*σ*_{0}) → *σ*/*σ*_{0} in the low shear rate limit as *η* → *η*_{0} in this limit. A specific example of Eq. (15) is the Eyring formula,^{8,9}

where *F*(*σ*/*σ*_{0}) = sinh(*σ*/*σ*_{0}). In this case, *σ*_{0} = *k*_{B}*T*/2$v0$, where $v0$ is a characteristic volume traced out by the flowing “unit” (again not well defined on a molecular scale) during a cage of the molecule reorganisation event or “jump.” By rearrangement of Eq. (16),

where *τ*_{0} = *η*_{0}/*σ*_{0} is a characteristic relaxation time and $\sigma =\eta 0\gamma \u0307$ has been used. Note this is not the Maxwell relaxation time, *τ*_{M} = *η*_{0}/*G*_{∞}, where *G*_{∞} is the infinite frequency shear rigidity modulus. The Newtonian viscosity acts as a coefficient or reference state whose temperature dependence is often assumed to follow an activation energy or Arrhenius dependence although there are a number of alternative assumed flow mechanisms proposed in the literature which also give rise to the same inverse temperature dependence but do not invoke the concept of activated dynamics. Notably there is the so-called “shoving model,”^{39–43} in which stress relaxation is taken to be caused by elastic deformation and ensuing structural rearrangement of the cage of molecules around an arbitrary molecule. The Arrhenius activation energy in the shoving model is replaced by the product of *G*_{p}, the plateau modulus (which is slightly lower than the true infinite frequency modulus, *G*_{∞}), and Ω_{c}, a characteristic volume of order of the molecular volume or smaller.

What is not so well known is that Eyring’s shear thinning formula can also be applied to predict time dependent shear rate viscoelastic effects such as during shear-start-up,

when it is incorporated in Maxwell’s equation including time which leads to viscoelasticity,^{44}

which expresses the total strain rate as the sum of elastic and viscous parts, the first two terms on the left of Eq. (19), respectively. Note that the infinite frequency shear rigidity modulus, *G*_{∞}, is assumed to be independent of $\gamma \u0307$ in the Eyring model, which is usually not a bad approximation.^{35}

Substitution of Eq. (16) in Eq. (19) gives the fundamental equation used in our extension of the EM,

The analytic solution of *σ*(*t*) from Eq. (20) is presented in Eq. (A8) derived in Appendix. Figure 4 shows what might be called a “time dependent instantaneous viscosity,” $\eta (t)=\sigma (t)/\gamma \u0307$, normalized by the Newtonian value, *η*_{0}, for several shear rates, starting from an equilibrium (i.e., *σ* = 0) state. Figure 5 is as is shown in Fig. 4 except that the starting stress at *t* = 0 corresponds to $\gamma \u0307=1.5$ and the shear rate for *t* > 0 is 1.5 + *x* where *x* is the number given in the figure. This results in some of the curves representing the effect of a decrease in shear rate at *t* = 0. The figure demonstrates that Eq. (A8) is capable of predicting the stress response for a step-down as well as a step-up in shear rate, and is one of the main results of this work.

The Transient Time Correlation Function (TTCF) for shear stress, $CTr(t,\gamma \u0307)$, is^{16,47,48}

where ⟨⋯⟩ is a time correlation function, which is the product of the stress taken from an equilibrium (i.e., unsheared) system at *t* = 0 and one at a time *t*, *σ*(*t*), which has evolved at constant shear rate, $\gamma \u0307$, during that period (note the shear stress is positive for, $\gamma \u0307>0$).

The stress build up after the commencement of shear is more complicated if expressed in terms of an alternative stress relaxation formulation formalism used in the rheological literature, in which the SRF is usually expressed as a weighted sum of exponentials with different relaxation times. The current relaxation times are a function of the present physical state of the system as it evolves from an equilibrium to a sheared steady state. If we assume that the current shear stress relaxation function at time *t*′ is *G*_{∞}*ϕ*_{s}(*τ*_{s}(*t*′)), where *τ*_{s} is a characteristic shear stress relaxation time, then the time dependent stress from shear start-up at an arbitrary shear rate would be^{45}

which is more complicated than the TTCF approach given in Eq. (21). The two treatments in Eqs. (21) and (22) are both correct, but the two formulations are in terms of different time dependent functions.

Figure 6 shows $CTr(t,\gamma \u0307)$, normalized by *G*_{∞}, for several $\gamma \u0307$ values using Eq. (A10) for the Eyring model. The function decays more rapidly with time with increasing shear rate. Substitution of Eq. (A9) in the definition of *η*_{i} in Eq. (2) gives for the Eyring model [i.e., where Eq. (17) applies],

which is, in fact, a special case of the Carreau shear thinning formula (i.e., with the exponent of the denominator set to 1/2),^{49} sometimes used to represent the viscosity, *η*, itself. An advantage of the incremental viscosity in the present case is that the Carreau model (with *α* = 1/2) has a point of inflection at a shear rate of 0.71/*τ*_{0} whereas the Eyring model has a point of inflection at 1.035/*τ*_{0}. This should facilitate a more accurate extrapolation of the shear rate dependent viscosity to zero shear rate and hence to the Newtonian viscosity (note that *η* and *η*_{i} have the same value in the $\gamma \u0307\u21920$ limit). The relationship between *η* and *η*_{i} for other simple analytic forms for *η* is discussed further in the supplementary material.

Figure 7 presents a comparison of *η* and *η*_{i} as a function of the shear rate for the Eyring model using Eqs. (17) and (23), respectively. The incremental viscosity, *η*_{i}, is seen to decay more rapidly with shear rate than *η*, as would be expected from the defining equations in Eq. (3).

The incremental viscosity can be calculated if an infinitesimally small shear rate “kick,” $\gamma \u0307k$, is added at *t* = 0 to an equilibrium or steady state sheared system at shear rate $\gamma \u0307$. Now, *σ*_{k}(*t*) is the difference in the shear stress between the kicked and reference trajectories and therefore *σ*_{k}(0) = 0. The viscoelastic response associated with the shear-kick can be written in a form involving a shear stress relaxation function appropriate to the $\gamma \u0307$ state,

Figure 8 shows the Eyring model incremental viscosity stress relaxation function, *C*_{i}, given in Eq. (A11) for several shear rates, revealing that it decays more rapidly with time with increasing (reference state) shear rate. The slower relaxation processes are more affected by the shearing mechanism, which is the same qualitative trend as found in past NEMD studies.^{50} Comparison between Fig. 6 and Fig. 8 reveals that the $Ci(t,\gamma \u0307)$ are more separated than the $Ctr(t,\gamma \u0307)$ for the same $\gamma \u0307$, notably at short times.

See the supplementary material for discussion of other model analytic forms for *η* found in the literature and their derived *η*_{i}.

Section III deals with the extent to which the NEMD data and Eyring model comply with the Fluctuation Theorem.

## III. FLUCTUATIONS AND THE FLUCTUATION THEOREM

The average stress, *σ*_{t}, over a single time period from 0 to *t* is

Let us assume that the set of *σ*_{t} follows a shifted gaussian probability distribution,

where $\sigma \xaf$ is the simulation average stress taken over many segments and *σ*_{D}, is the standard deviation (SD) of the quantity, $[\sigma t\u2212\sigma \xaf]$, over that period ($>>$*t*). Then,

For an equilibrium system,^{14,51} the Einstein-Helfand formula for the viscosity ($i.e.,\u2009\sigma \xaf/\gamma \u0307$) is, for times *t*, greater than the Maxwell relaxation time, with $\sigma \xaf=0$,

where ⟨⋯⟩ indicates an average over a sequence of segments of length *t*. We assume here that it is still reasonably accurate for finite but not too large $\sigma \xaf$. Equation (28) substituted in Eq. (27) (assuming the shifted gaussian stress PDF is applicable at a finite shear rate) gives the Fluctuation Theorem (FT) expression,

Despite making a number of assumptions, the above short approximate derivation of the FT for shear does lead one to expect that a (shifted) gaussian distribution of time-averaged shear stresses should still be compatible to a reasonable approximation with the FT. Note however that the FT does not rely on the PDF being gaussian. The FT is not able to provide the analytic form of $\eta (\gamma \u0307)$ but any model system that gives rise to a particular rheological constitutive equation should comply with the FT. Although Eq. (27) is only formally exact in the equilibrium ($i.e.,\u2009\gamma \u0307\u21920$) limit, it is a reasonable approximation for not too high shear rates, especially for long segment averaging times.^{52,53} Deviations from the gaussian form, however, have important consequences which, for example, limit the use of the Green-Kubo expression for the viscosity to small values of the shear rate.^{52}

### A. NEMD results and the Fluctuation theorem

A spatially local stress can be defined, and the effect of the stress sampling volume on the PDF for finite shear rates is investigated here, to explore what the analytic form of this stress distribution actually is. Basing this definition on Eq. (12), the contribution from molecule *i* to the virial expression for the shear stress of the whole periodic system is

where the average number density is *ρ* = *N*/*V*. The virial expression for the stress of the whole system is $\sigma =\u2211i=1N\sigma i/N$ and for a subvolume Ω in the simulation cell, containing on average *N*_{Ω} molecules,

which is the average of *σ*_{i} for those molecules inside a smaller volume Ω of a chosen shape. This Atomic Virial (AV)^{54} definition is just one of many possible prescriptions for the local stress. It has the advantage that it is physically reasonable, being based on those molecules located in the subvolume, and it goes over to the usual virial expression when Ω → *V*, the volume of the simulation cell. The AV method is suitable for spatially homogeneous systems but gives rise to unphysical oscillatory stress fluctuations near walls, in which case the exact Method of Planes (MOP) method^{55} and Volume Averaging (VA)^{56} methods should be used.

The logarithm of the probability distribution functions (PDF) of *σ*_{Ω} for equilibrium and sheared liquids is presented in Fig. 9 using a spherical subvolume, Ω (see Refs. 54 and 57 for further details). The average number of molecules in the subvolume is indicated in the figure. The $\gamma \u0307=0$ distributions increasingly depart from a gaussian form as the sampling volume or *N*_{Ω} decreases. For not too small *N*_{Ω}, the main effect of the shear rate is to shift the gaussian along the abscissa by the average shear stress, $\sigma \xaf$, value. For smaller *N*_{Ω}, the distribution adopts a noticeable skewed peaked shape with prominent tails in the high stress wings. These trends are the same for a cubically shaped subvolume. Figure 10 shows the PDFs as a function of the shear rate for spherical and cubic subvolumes of volume, 4*πR*^{3}/3 and 8*R*^{3}, respectively, where *R* = 2.41. The gaussian shape is reasonably well maintained but the standard deviation increases with $\gamma \u0307$. This partly vindicates an early phenomenological molecular theory of viscosity for simple liquids which made this assumption.^{58}

The shape of the PDF can be characterized using the second to fourth moments of the shear stress distribution about the mean, *σ*_{2}, *σ*_{3}, and *σ*_{4}.^{59} The standard deviation (SD), $\sigma 2$, skew (SK), $=\sigma 3/\sigma 23/2$, and kurtosis (KU), $X4=\sigma 4/\sigma 22\u22123$, measure the breadth, degree of slant, and extent of the tails, respectively, of the PDF. Figure 11 shows that these quantities increase with decreasing Ω, which is consistent with the PDF shapes seen in Fig. 9. The SD, SK, and KU increase in magnitude with diminishing subvolume, to a good approximation as ∼Ω^{−ν}, where *ν* is an exponent. There is a gradual departure from a gaussian form as the size of Ω decreases. Searles and Evans showed that the instantaneous and time averaged shear stress for the whole of a sheared liquid is close to being gaussian.^{52} The trends in Fig. 11 are consistent with the higher moments being larger than those of a gaussian.^{60–64} The effect of shear is also to skew the PDF to the right. The figure also shows that the data for sphere and cubes fall on the same lines, so the magnitude of the subvolume, Ω, rather than its shape (within reasonable limits) is the key governing factor for this local stress model.

### B. The Eyring model and the Fluctuation theorem

The EM assumes that, on the one hand, there is only one shear stress and one shear rate in the system, while on the other hand, the flow process involves molecular “jumps” which would imply strong local fluctuations in both quantities. The two scenarios sit rather uncomfortably within the same model. This issue is considered in this subsection and partly resolved. A modification of the EM is derived in which the effects of local shear stress fluctuations are incorporated through the introduction of a gaussian shear stress PDF.

Let the time and spatial average value of the local shear stress, *σ*_{t}, be $\sigma \xaf$. The EM can be generalized by assuming that there is an assembly of such local states at any given time characterized by *σ*_{t}, each obeying

where *A*_{u} = *σ*_{0}/*η*_{0}. The *σ*_{t} distribution is taken to be gaussian from Eq. (26). The system is deemed to explore a distribution of *P*(*σ*_{t}) and hopping rates, which is in fact the Ree-Eyring extension to EM but with the average flow a sum of stressed units, rather than the other way around.^{65,66} Previous shear stress relaxation models^{67–69} have assumed such processes, which are also consistent with stretched exponential stress relaxation.^{70} There is then a distribution of shear rates at any one time, whose average value over time and space is

The latter step in Eq. (33) was performed by integrating the two exponential terms in the sinh function separately and using $\u222b\u2212\u221e\u221e\u2061exp(\xb1ax)exp(\u2212b(x\u2212c)2)dx=\pi /bexp([a2/4b]\xb1ac)$.^{71} In the zero shear rate limit, from Eq. (33),

to give $Au=exp(\u2212\sigma D2/2\sigma 02)\sigma 0/\eta 0$, which when combined with Eq. (33) leads to

This is the same analytic form as the original Eyring relation of Eq. (32) with the single stress, *σ*_{t}, replaced by the system’s average value, $\sigma \xaf$. Therefore, even though there is a stress PDF which is much broader than the (Eyring) delta function distribution, the analytic form of the final formula relating average shear stress and shear rate is the same. The fact that the Eyring analytic form is consistent with a realistic distribution of stresses on the molecular scale might in part explain its success in reproducing finite shear rate experimental data, as *σ*_{0} is almost always treated as a parameter to be fitted to experimental or simulation data.

### C. Microstructural origins of shear thinning in NEMD and the Eyring model

The above treatment does not indicate how large these flowing units are. Insights into this can be gained for the present LJ system by considering a spatially resolved shear stress correlation function,^{72–74}

which measures correlations between the configurational stress experienced by two particles *i* and *j* a distance *r* apart. The kinetic component of the stress is not included in the above definition, but for liquid-like densities, and not too high *T*, it is relatively small compared to the interaction part of the stress. Unlike *g*(*r*), it is defined here so that it is not normalized by 4*πr*^{2}Δ*r* and as a quantity is more compatible with the number of particles in a shell of thickness Δ*r* at *r*. Figure 12 shows *P*_{xy}(*r*) for a dense equilibrium fluid state (the simulation details are given in the figure caption). Note that if *P*_{xy}(*r*) = 0, there is no shear stress correlation between particles *i* and *k* at that separation. The figure shows that stress correlations reflected in the oscillations extend to at least the third coordination shell as found also in Refs. 72 and 73, which contains of order 100 particles. The number of particles within a sphere of radius *r*, or *n*(*r*), is also given in the figure. Such long range correlations in the local shear stress might have been anticipated from the work of Stassen and Steele who decomposed the stress time correlation function used in the Green-Kubo formula for the Newtonian viscosity into its two-, three-, and four-body components.^{75} They were all shown to be important and to decay slowly at long times.^{75} This behavior also explains why there is a strong system size dependence of the shear viscosity of high density fluid states, especially near the triple point.^{73} Figure 12 suggests that for simple monatomic liquids, the number of molecules involved in the flow process is rather large (∼100), which is probably much larger than would be the case for large organic molecules. The relative success of the Eyring model for large molecules is possibly in part attributable to a more localized stress field than for simple monotomic liquids.

Simple molecular fluids (e.g., Lennard-Jones) are one type of system where the Eyring equation has had mixed success in predicting the shear thinning curve.^{6,76} This is illustrated using the present data in Fig. 13, which shows the shear thinning curve for a LJ system close to the triple point (note the log-log scale). The Eyring model does not match the simulation data for large shear rates, where a power law is better. This is possibly because the shear rates are relatively high (∼10^{12} s^{−1} for liquids composed of small molecules such as argon) for LJ molecules and long range structural changes induced by the shearing action are more readily achieved due to their spherical shape, which is expressed in its most extreme form at very high shear rates by the appearance of the so-called “string” phase.^{77} The possibility of long-range reordering due to shear was not envisaged in the construction of the Eyring model, only local small shear deformations of the structure. In addition, for a LJ liquid or real liquid composed of small molecules, a significant proportion of the area under the shear stress relaxation function (and hence viscosity) at short times is due to kinetic or inertial stress decay, whose time expansion is in even powers of time [unlike the Eyring model stress relaxation decay which is based on exponentials; see $\sigma (t,\gamma \u0307)$ in Eq. (A8) and derived functions].

Another limitation of the Eyring model is that it assumes that for large shear rates, the stress continues to increase with shear rate, whereas in some studies (e.g., for glasses^{78}), a limiting stress value has been found.

### D. Fictive Eyring

Another limitation of the Maxwell-Eyring model of viscoelasticity is that it does not predict shear stress overshoot (where the shear stress exceeds at short time its long time value) at high shear rates and soon after the start-up of shearing, an effect observed in polymer solutions, melts,^{80} metallic glasses,^{79} and even in simple liquids.^{50} This absence is because the instantaneous relaxation time soon after the application of a step in the shear rate is calculated from the steady state shear formula [see Eq. (20)], which is questionable as the system may not have had time to come to a structural and rheological steady state at that shear rate during this transient period. A real system takes time to evolve its structure (and hence stress) to that of the final steady state. Therefore at short times, the EM underestimates the structural relaxation time. The real system’s characteristic structural and property values would lag behind those predicted using the Maxwell-Eyring viscoelastic model, possibly resulting in a shear stress overshoot.

In the field of glass formation, a so-called “fictive” temperature, *T*_{f}, has been introduced to quantify such time lag effects. The fictive temperature is that of an equilibrium liquid which has the same structure and thermodynamic properties as a supercooled state. The fictive temperature is higher than or equal to the actual temperature of the system, i.e., *T*_{f} ≥ *T* during rapid cooling.^{81} This takes account of the lag in structural evolution as the system is forced out of equilibrium.

A similar construction can be adopted for the shear stress evolution as introduced in Ref. 45 to represent the physical state of a lubricant fluid element as it passes through the contact zone in high pressure elastohydrodynamic lubrication. The stress relaxation function in that study was represented by a stretched exponential in time, which can be written as a weighted sum of exponentials with a distribution of relaxation times. In the present simpler analytic model, there is only one effective relaxation process representing the system at any given time after shear-start-up [see Eq. (19)]. This time lag effect can be introduced empirically in the Eyring model by the inclusion of a time-delay function, $Fi(t,\gamma \u0307)$, which has the limits of 0 and 1 at *t* = 0 and *t* = *∞*, respectively. From Eq. (19),

The performance of one form chosen for $Fi(t,\gamma \u0307)=1\u2212exp(\u2212\u03f5\gamma \u0307t)$, where *ϵ* is a positive constant, is seen in Fig. 14. The figure shows the case where shear start-up from an equilibrium liquid is initiated with $\gamma \u0307=1$ and *ϵ* = 4 and using the same Eyring parameters as for the previous figures. The functional form of $Fi(t,\gamma \u0307)$ and value of *ϵ* are physically reasonable as the peak in the overshoot stress profile is typically in the strain ($i.e.,\u2009\u222b0t\gamma \u0307(x)dx$) range between 0.1 and 1.0. If *F*_{i} = 1, we have the original Eyring model, and for the adopted time-delay functional form, this case might be called the *fictive Eyring* model. Although there is no specific fictive variable in Eq. (37), its effects are encapsulated in the additional semi-empirical term, *F*_{i}. Figure 14 shows that this model exhibits a peak in the stress at short times, unlike the original Eyring model itself (i.e., where $Fi(t,\gamma \u0307)=1$) also shown in the figure, which *always* increases monotonically with time after shear start-up. An alternative approach to generate a time delay effect would be to relate the instantaneous relaxation time to an order parameter which evolves more slowly than, in this case, the Eyring shear stress. For temperature changes in the glassy state, a fictive temperature has been used to implement this time delay.^{82,83} In Ref. 83 this was interpreted in terms of a so-called material time which is an effective time which is moving slower than the actual time. For a pressure jump, a fictive free volume has been used.^{22} In the present case, this order parameter would be evolved parallel in time with the shear stress itself. Further experimental and NEMD studies would be needed to identify a suitable candidate for this property.

## IV. CONCLUSIONS

Rheometers measure shear stress against shear strain or strain rate, $i.e.,\u2009\sigma (\gamma \u0307)$, and the invocation of a “shear viscosity” may not be a necessary or desirable exercise as it is certainly not uniquely defined apart from in the Newtonian regime where the shear rates are very small on a molecular time scale. The usual definition of the non-Newtonian viscosity is the ratio of the shear stress divided by the shear rate ($i.e.,\u2009\eta (\gamma \u0307)=\sigma (\gamma \u0307)/\gamma \u0307$). We propose an alternative definition of a shear rate dependent viscosity, the derivative of shear stress against strain rate, which we call the “incremental viscosity” (IV) or $\eta i(\gamma \u0307)$. The IV is a better representation of the physical state of the system at a particular shear rate than $\eta (\gamma \u0307)$. Also *η*_{i} converges to the Newtonian viscosity as the shear rate tends to zero, and because the IV is compatible with Boltzmann’s superposition principle, it is the natural extension of the Newtonian viscosity definition to systems where the shear rate is arbitrarily large. The IV is calculated here by NEMD using a time dependent shear rate or “shear-kick” scheme and analytically for the Eyring model where it was found to be a special case of the Carreau formula which has often been used to represent liquid $\eta (\gamma \u0307)$. The associated stress relaxation function for the shear kick is also computed and derived within the NEMD and Eyring models.

The Eyring model formula is shown to apply also when a gaussian distribution of stressed regions in space is assumed, which is approximately compatible with the fluctuation theorem. This may in part explain the Eyring model’s success in reproducing experimental data for a wide range of liquid classes. A simple semi-empirical modification of the Eyring model for the time dependent stress after shear start-up, to include the effects of time delay in the system’s response, is also shown to lead to shear stress overshoot behavior which is often seen in glassy and polymer solution systems.^{80,81}

The concept of the incremental viscosity unifies the rheological description of both Newtonian and non-Newtonian regimes, as the viscosity in both cases is defined as the derivative of the stress with respect to shear rate. The Newtonian viscosity has to be defined that way as any shear rate dependence of the ratio of stress to strain rate (i.e., not a linear relationship) means that one is not in the Newtonian shear rate regime. In practice, the incremental viscosity could be obtained experimentally by the direct application of a small shear rate pulse or an oscillatory strain rate to an already sheared system, as in the NEMD calculations reported here. Alternatively, the incremental viscosity could be computed by fitting experimental or simulation shear stress *vs.* shear rate data to a function and then analytically differentiating this with respect to shear rate.

We would like to add that the usual definition for the viscosity ($i.e.,\u2009\sigma /\gamma \u0307$) has advantages, for example, in being directly measurable by experiment or NEMD simulation, and the stress of the system is readily obtained from it. The incremental viscosity has complementary advantages, however, in probing more accurately the physical state of the system at a given shear rate and in being compatible with viscoelastic constitutive equations based on Boltzmann’s superposition principle.

## SUPPLEMENTARY MATERIAL

See supplementary material for further discussion of the relationship between *η* and *η*_{i} for other simple analytic forms for *η* found in the literature.

## ACKNOWLEDGMENTS

D.D. would like to acknowledge the support received from the EPSRC under the Established Career Fellowship Grant No. EP/N025954/1. All data can be made available by emailing the authors of the paper or tribology@imperial.ac.uk.

### APPENDIX: DERIVATION OF EYRING VISCOELASTIC FORMULAS

The analytic solution of Eq. (20) for *σ*(*t*) is given here.

Let *t* = 0 be the time an equilibrium liquid is subjected to a constant shear rate. An analytic solution of Eq. (20) was derived in Ref. 45 via a Riccatti equation to give

where

An alternative more general solution for *σ*(*t*) resulting from the following strain rate history,

is derived here. This represents the transformation of the system with time from a steady stressed state with shear rate, $\gamma \u03071$, at *t* = 0^{−} to eventually a new steady state at $\gamma \u03072$, of arbitrary relative magnitude and sign. This is derived using the substitution, *Y*(*t*) = exp(*σ*(*t*)/*σ*_{0}) in Eq. (20), which after rearrangement gives

where *A*_{1} = −*G*_{∞}/2*η*_{0}, $B1=\gamma \u0307G\u221e/\sigma 0$, and *C*_{1} = *G*_{∞}/2*η*_{0} are time independent parameters. Equation (A4) can be integrated by separation of variables to give^{46}

where

and *G* is the constant of integration. Therefore, if *H* = exp(*SG*),

where *Y*(0) = exp[*σ*(0)/*σ*_{0}] is the value of *Y* at *t* = 0, the starting point being a steady state system with stress, *σ*(0) = *σ*_{1}. From Eqs. (A5) and (A7),

which is the time dependent stress during the transition between the two steady states 1 and 2. Equation (A8) has the long time limit,

which reduces to the Eyring formula of Eq. (17) on substitution of the definitions of the variables *A*_{1}, *B*_{1}, and *S* in Eq. (A9). Note in Eq. (A9) that *A*_{1} is negative, while *B*_{1} and *S* are >0.

The Transient Time Correlation Function (TTCF) for shear stress $CTr(t,\gamma \u0307)$ defined in Eq. (21) is within the Eyring model approximation, $CTr(t)=(d\sigma (t)/dt)/\gamma \u0307$ [in analogy to Eq. (11)], which from Eq. (A8) is

where *σ*(0) = 0. In the Eyring model, *G*_{∞} is independent of shear rate, which is why $CTr(t,\gamma \u0307)$ is conveniently defined in Eq. (A10) by normalizing it by *G*_{∞}. The reason the shear stress TTCF in the Eyring model is relatively straightforward to derive is that this time dependent function depends only on the *instantaneous* stress at time *t*′ where 0 ≤ *t*′ ≤ *t*, i.e., without any time delay.^{29} It is essential to know *σ*(*t*) by some other means such as from Eq. (A8) or directly by NEMD in order to calculate the TTCF.

Substituting Eq. (A8) into Eq. (24), an analytic expression for the incremental viscosity relaxation function, *C*_{i}(*t*) within the Eyring model is

Just as for *C*_{Tr} in Eq. (A10) and *C*_{Tr}, the relaxation function is normalized by *G*_{∞}, which is assumed to be independent of shear rate. The difference between the Eyring formula for the TTCF in Eq. (A10) and *C*_{i} in Eq. (A11) is that for *C*_{i}, the relevant stress is the difference between that of the $\gamma \u0307+\gamma \u0307k$ and $\gamma \u0307$ trajectories, and the normalizing factor in Eq. (A11) is $\gamma \u0307k$ rather than $\gamma \u0307$ in Eq. (A10). When $\gamma \u0307\u21920$, the two functions are identical.